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I.  INTRODUCTION 


This  final  report  covers  research  carried  out  in  conjunction  with  the  Plasma  Physics 
Division  of  the  Naval  Research  Laboratory  for  the  p^yformance  period  10  December  1984 
to  9  December  1985  for  Contract  N00014-86-C-2069.  The  research  involved  the  theoretical 
and  numerical  analysis  of  the  physics  of  free  electron  lasers  using  relativistic  particle  beams. 
The  major  emphasis  of  the  research  was  to  obtain  design  criteria  for  the  development  of 
a  two-stage  FEL  oscillator  operating  in  the  trapped  particle  mode.  The  present  work  was 
centered  on  the  development  of  a  fully  relativistic,  nonlinear  analysis  of  the  spatial  and 
temporal  evolution  of  multiple  modes  within  a  free  electron  laser  oscillator  and  of  a  large 
amplitude,  nonlinearlv  saturated  state  characteristic  of  trapped  particle  mode  operation. 


The  equations  solved  are  the  Maxwell  equations  of  electrodynamics  coupled  with  the 
collisionless  Boltzmann  equation  that  describes  collisionless  particles  under  influence  of  the 
electromagnetic  fields.  The  electromagnetic  fields  include  the  radiation  fields  from  the  FEL 
and  the  self-electric  fields  from  the  longitudinal  potential  due  to  the  space  charge,  i.e.,  the 
dominant  component  of  the  interparticle  Coulomb  forces.  The  particle  dynamics  transverse 
to  the  magnetic  axis  are  included,  but  gradients  in  the  radiation  fields  are  ignored.  The 
electron  beam  equilibrium  is  assumed  to  be  spatially  uniform  and  temporally  stationary. 
Justification  and  probable  impact  of  further  approximations  are  discussed  in  the  technical 
section  of  this  report.  The  approximations  employed  are  consistent  with  the  purpose  of 
obtaining  experimentally  implementable  design  criteria  for  the  FEL  oscillator. 

Both  analytical  and  numerical  analyses  were  performed.  The  numerical  simulations 
are  in  qualitative  and  quantitative  agreement  with  the  analytical  theories.  For  the  FEL’s  of 
interest,  the  theory  exhibits  the  same  scaling  as  is  obtained  for  an  FEL  amplifier  operating 
in  the  low  gain  regime  and  the  ultrahigh  gain  regime.  The  intermediate  case,  the  moderate 
gain  regime,  is  directly  applicable  to  the  experimental  parameter^  of  the  FEL  at  NRL.  For 
example,  for  a  beam  energy  of  500  keV,  a  current  of  100  A,  radius  of  0.64  cm  and  a  wiggler 
length  of  4.0  cm  with  wiggler  field  strength  of  615  G,  the  theoretical  expression  for  the 
threshold  reflection  coefficient  is  0.64.  The  experimentally  measured  value  is  0.65,  a  very 


9 


satisfactory  agreement.  This  and  other  examples  may  be  found  in  the  NRL  Memorandum 
Report  5G79  contained  in  Appendix  B.  This  NRL  memorandum  report  has  also  been 
published  in  Nuclear  Instruments  and  Methods  in  Physics  Research  A250,  159-167  (19S6). 
Appendix  A  contains  the  FORTRAN  listing  of  the  computer  codes  that  were  written  for 
use  in  the  research. 
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II.  TECHNICAL  DISCUSSION 


We  have  conducted  an  analytical  and  numerical  analysis  of  the  field  evolution  in  a 
high  gain  free  electron  laser  (FEL)  that  is  operating  in  the  oscillator  configuration.  This 
high  gain  oscillator  provides  the  pump  field  for  the  second  stage  of  a  two  stage  free  electron 
laser  that  is  operating  in  the  trapped  particle  regime. 

The  analysis  contains  the  following  attributes.  The  electron  beam  is  assumed  to  be 
spatially  uniform  and  temporally  stationary.  The  magnetostatic  wiggler  field  is  helically 
symmetric  and  generates  circularly  polarized  radiation  with  the  same  polarization  sense 
as  that  of  the  electron  beam  motion  in  the  magnetic  field.  The  motion  of  the  electrons 
in  the  combined  radiation  and  magnetic  fields  produces  ponderomotive  bunching  which 
is  inhibited  by  the  self  consistently  obtained  space  charge  fields  which  arise  due  to  the 
particle  bunching. 

In  order  to  obtain  a  detailed  understanding  of  the  distinct  physical  processes  which 
contribute  to  the  operation  of  the  FEL  in  this  configuration,  several  specific  computer 
codes  were  developed  and  are  summarized  in  the  following  discussion. 

The  first  of  these  is  entitled  MULTI.FOR  and  enables  one  to  monitor  the  simultaneous 
evolution  of  many  longitudinal  modes  of  the  radiation  field.  These  fields  are  given  by 

Ar(z,  t)  =  o,n{t)sin{knz)exp{iunt)e-  +  c.c.  ( 1) 

n 

and 

4>{z,  t)  =  <t>insin[(kn  +  kw)z  -  u>nt ]  +  <f>2ncos[(kn  +  kw)z  -  unt\  (2) 

n 

where  A n(z,t)  is  the  vector  potential  for  the  circularly  polarized  radiation,  <f>(z,t)  is  the 
self-consistent  longitudinal  space  charge  potential,  c.c.  is  the  complex  conjugate  of  the  first 
term  on  the  right  of  Equation  (1)  and  kw  is  the  wiggler  wave  number.  Computationally, 
the  decomposition  of  the  fields  can  contain  at  most  fifty  modes.  The  operation  of  this  code 
yields  information  on  the  growth  properties  and  saturated  values  for  the  radiation  field. 
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For  the  sake  of  completeness,  computer  codes  were  developed  to  evaluate  the  the¬ 
oretical  expressions  for  the  growth  rate  for  comparison  to  the  simulation  results.  Two 
computer  codes  were  developed  for  this  purpose  and  they  are  entitled  PARTEMP.FOR 
and  MAXWTEM.FOR.  In  addition  to  being  capable  of  evaluating  the  growth  rate  for  the 
cold  beam  case,  these  codes  have  the  capacity  to  evaluate  the  growth  rate  when  the  elec¬ 
tron  beam  possesses  a  spread  in  parallel  energy.  In  the  code  PARTEMP.FOR,  the  parallel 
temperature  spread  is  modeled  with  a  Lorentzian  distribution  in  axial  momentum  and  the 
associated  dispersion  relation  is  evaluated  by  quadratic  interpolation  methods. 

The  code  MAXWTEM.FOR  models  the  axial  spread  with  a  Maxwellion  distribution. 
The  associated  dispersion  relation  is  expressed  in  terms  of  the  Fied-Conte  plasma  dispersion 
function  and  the  growth  rates  are  again  evaluated  by  quadratic  interpolation  methods.  By 
defining  a  characteristic  axial  spread,  A pz,  which  encompasses  70%  of  the  beam  electrons, 
one  notes  that  irrespective  of  the  detailed  structure  of  the  axial  distribution  one  obtains 
the  same  growth  rate.  Now,  comparing  this  growth  rate  to  a  simulation  case,  as  shown  in 
Figure  1.  One  notes  that  the  horizontal  line  representing  the  theoretical  growth  rate  for 
the  mode  under  consideration  is  in  excellent  agreement  with  the  simulation  results  for  the 
linear  gain  regime  for  which  the  comparison  is  appropriate. 

For  significant  reductions  in  the  electron  beam  current  or  wiggler  field  strength  the 
growth  rate  for  the  generated  radiation  is  also  reduced.  These  conditions  require  calcula¬ 
tions  with  increased  accuracy,  and  for  this  reason  a  double  precision  version  of  the  field 
evolution  code  was  developed.  This  code  is  entitled  SPRADB.FOR  and  returns  all  the 
capabilities  of  the  single  precision  code  MULTI. FOR. 

The  increased  storage  requirements  for  the  double  precision  variables  in  addition  to 
the  increase  cpu  time  for  operations  on  these  variables  made  long  runs  or  large  numbers 
of  particles  prohibative  for  code  operation  on  the  VAX  780  or  785.  For  these  reasons  a 
Cray  FORTRAN  code  was  developed  and  entitled  GRAYWIG4.FOR.  This  code  made  use 
of  the  Cray  processing  advantages  while  retaining  the  attributes  and  capabilities  of  the 
code  developed  on  the  VAX.  Further  optimization  of  the  code  data  handling  capabilities 
were  incorporated  in  the  Cray  codes  entitled  CMPLXAMP.FOR  and  RE3TART.FOR. 


The  code  COMPLXAMP.FOR  made  use  of  the  complex  amplitude  representation  of  the 
radiation  fields  to  evolve  the  system  of  equations.  This  method  proved  to  be  faster  and 
more  stable  than  its  Cray  predecessor.  The  final  code,  RESTART. FOR,  was  developed  to 
take  previous  results  from  CMPLXAMP.FOR  as  input  in  order  to  evolve  the  system  for 
longer  times. 

The  application  of  these  computational  tools  in  the  analysis  of  experimental  results 
is  presented  in  the  attached  publication  in  Appendix  B.  The  quantitative  and  qualitative 
comparisons  are  excellent. 
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File  ^VC$DRBl : [ MARABLE ] CMPLXAMP . CFT ; 1  (4034,5,0),  last  revised  on  21-MAR-1988  14 
: 0 5 ,  Is  a  34  block  sequential  file  owned  by  UIC  [MARABLE).  The  records  are  vari 
able  length  with  implied  (CR)  carriage  control.  The  longest  record  is  72  bytes. 


Job  CMPLXAMP  (1995)  queued  to  LN03_QUE  on  21-MAR-1988  14:08  by  user  MARABLE, 
UIC  (MARABLE],  under  account  4790  at  priority  100,  started  on  printer  LTA8 :  on 
21-MAR-1988  14:08  from  queue  VCLN03B. 
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CODE  TO  EVALUATE  TEMPORAL  EVOLUTION  OF  THE  SPRECTRA 
OF  UNSTABLE  MODES  IN  A  HELICAL  WIGGLER  FREE  ELECTRON  LASER 
DELETION  OF  FIRST  TRANSIT  TIME 

FIELD  AND  PARTICLE  EQUATIONS  ARE  EVOLVED  BY  ADAMS-BASHFORTH 
METHOD  WITH  INITALI ZATION  BY  RUNGE-KUTTA  METHOD 
REFORMULATION  OF  THE  PARTICLE  PHASE  3/13 
CONVERSION  TO  CRAY  FORTRAN 

INCLUSION  OF  FREQUENCY  SHIFT  ERROR  CHECK  IN  ADAMS-BASHFORTH 
EQUATION  SOLVER 

REPLACE  EXPRESSION  FOR  THE  DERIVATIVES  WITH  THE  FUNCTIONAL 

EVALUATION  OF  THE  DIFFERENTIAL  EQUATION 

COMPLEX  AMPLITUDE  FORMULATION  OF  FIELD  EVOLUTION 

REAL  BETAZO , BETAO , KPOD( 20 ) ,OMEGA(20) , BETAZ ,  GAM 

REAL  PSI(20,3000, 5) ,U0( 3000,5) 

REAL  TEMP( 3000) ,TEMP1( 3000) 

REAL  TIME( 5000) , PLOTA( 5000,20) , GROWTH ( 5000,20) 

REAL  FREQ( 5000,20) ,EWAV{5000) ,GROE(5000) 

REAL  PHI 1(20,5) ,PHI2(20,5) , PHI IT (2 0,4) ,PHl2T(20,4) ,KWIGL 
, PLAI ( 5000,20) , AMAG (20) 

REAL  KWIGR , NU, NUR , NUI , FILL, PLEE( 5000,3) ,PLAR( 5000,20) 

REAL  PSIT(20, 3000, 4) ,U0T( 3000,4) 

REAL  K31( 3000) ,K32( 3000) ,K3 3 ( 3000) ,K34( 3000) 

REAL  K41 (  3000 ) ,K42(3000) , K4 3 (  3000 ) , K44 (  3000 ) 

COMPLEX  ATEMP( 20 ) , A(20 , 5) ,APP( 20, 5) ,CORA( 20) ,AT( 20, 4 ) 

COMPLEX  APT (20, 4) , Kl 1(20), K1 2(20), K13( 20 ),K14(20) 

COMPLEX  FI 1(20), Fl 2(20), Fl 3(20), FI 4(20), Fl 5(20) 

COMPLEX  ALPHAl , CAPP ( 20 ) 

INTEGER  F , MM ,JP,J,J1,K, MAX I T 
COMMON/BLKl/BETAZO , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK2/PSIT , U0T , PHI IT , PHI 2T , AT , APT 

COMMON/BLK 3/TIME , NPART , N , RI SE , TAU , BETA1 , BETA2 , F , ALPHAl 
COMMON/BLK4/  ALPHA2 , ALPHA3 , NMODE 
COMMON/BLK 5/PS I , U0 , PHI 1 , PHI 2 , A , APP 
PARAMETER  ( PI=3 . 1415926535 ) 

PS00(J1,J)  =  -OMEGA( Jl) * FLOAT (J-1)*TAU 
TRISE(N)  -  1.  -RISE* ( EXP ( (— TIME(N)+1. ) /RISE )  - 
EXP ( -TIME ( N ) /RISE )  ) 

OPEN ( UNI T=1 , FI LE= ' FT01 ' , FORM= ' UNFORMATTED ' , STATUS= ' NEW ' ) 

OPEN( UNIT=2 , FILE=' FT02 ' , FORM= ' UNFORMATTED ' , STATUS= ' NEW ' ) 

FORMAT (  '  THE  PRED.-CORR.  METHD  FAILED  TO  CONVERGE  ON  STEP',2X, 

14,'  AFTER  '  , 1 4 , 2X ,  '  INTERATIONS '  ) 

WRITE ( 6 , 10 ) 

FORMAT (  'INPUT  NO.  OF  PART., NO.  OF  ITERAT ., GAM , BETAWIG , KWIGR 
, BUDKER , NWIG , EPS , PHASE , RISE , NPLUS , NMODE , NSEP , REF , F , FILL 
, MAXIT, ERROR, ERROR2 ' ) 

READ ( 5 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , NWIG , EPS , PHASE 
, RI SE , NPLUS , NMODE , NSEP , REF , F , FI LL , MAXI T , ERROR , ERROR2 
FORMAT ( '  INPUT  DATA : NPART , NTIMES , GAM0 , BETAW, KWIGR , NU , NWIG , 

EPS , PHASE , RISE , NPLUS , NMODE , NSEP , REF , F , FI LL , MAXIT , ERROR , ERROR2 ' ) 
WRITE( 6,20) 

WRITE( 6 , * ) NPART, NTIMES ,GAM0 , BETAW , KWI GR , NU , NWIG , EPS , PHASE , 

RI SE, NPLUS, NMODE, NSEP, REF, F, FILL, MAXIT, ERROR, ERROR2 
KWIGL  =  2 . * FLOAT ( NWIG) *PI 
BETAO  =  SQRT (1.0  -  1 . 0/( GAM0 *GAM0 )  ) 

BETAZO  -  SQRT( BETAO*BETA0  -  BETAW*BETAW ) 

NOPT  =  NINT( 2 . *FLOAT( NWIG ) * BETAZ 0/( 1 . -BETAZO ) )  +NPLUS 
OMEGA ( 1 )  =  ( FLOAT( NOPT ) *PI ) /BETAZO 
KPOD ( 1 )  =  KWIGL  +  BETAZ 0*OMEGA( 1 ) 

DO  50  J=2 , ( NMODE- 1 ) /2  +  1 

OMEGA ( J )  =  FLOAT ( NOPT+ (J-1)*NSEP) *PI/BETAZ0 
KPOD(J)  =  KWIGL  +BETAZ0  * OMEGA ( J ) 

OMEGA( NMODE+2-J )  =FLOAT( NOPT- ( J-l ) *NSEP ) *PI/BETAZ0 
KPOD( NMODE+2- J )  -  KWIGL  +  BETAZ0*OMEGA( NMODE+2-J ) 

CONTINUE 

BETA1  «  2 . *FILL*NU* KWIGL* *2*BETA0* BETAW/ ( KWIGR**2*BETAZ0**3 ) 

NUR  =  ( 1 . -REF )/BETAZ0 


c*** 

C*** 

c*** 


60 

C*  *  * 

80 

C*** 

90 

85 


1100 

1200 


NUI  =  -4 . *FILL*NU*KWIGL*  *2* ( 1 . -BETAW*  *2/2 .  )/( GAM0*BETAZ0  * 

1  KWIGR**2*BETAZO*OMEGA( 1 ) ) 

BETA2  =  8 . *NU*BETA0*KWIGL**2/( BETAZ0*KWIGR**2 ) 

BETA2  =  0. 

ALPHAl  -  -(NUR  +  CMPLX( 0  .  ,  1 .  ) *NUI ) 

ALPHA2  =  KPOD( 1 )/( BETAZO *BETAZ0 ) 

ALPHA3  =  . 25*KPOD( 1 ) *BETAW*GAM0/BETAZ0**2 
INITIALIZE  PHASE  AND  AMPLITUDE 
TIME ( 1  )  =  1. 

GROE(l)  =  0. 

TAU  =  1 ./FLOAT (NPART  -1) 

TAUT  =  FLOAT ( F ) *TAU 

INITIALIZE  EACH  MODE  UNDER  CONSIDERATION  AND  FIND 
THE  WAVE  ENERGY  DENSITY  IN  THE  INITIAL  SPECTURM 
EWAV0  =  0. 

DO  60  Jl-1, NMODE 
FREQ ( 1 , Jl )  =  0. 

GROWTH ( 1 , Jl )  -  0. 

PLOTA ( 1 , Jl )  =  EPS 
AMAG(Jl)  =  EPS 

A( Jl , 5 )  =  -CMPLX( 0  .  ,  1  .  ) *EPS*CEXP (  CMPLX ( 0 . , 1 . ) *PHASE  ) 

AT ( Jl , 1 )  =  A{ Jl , 5 ) 

EWAV0  =  EWAV0  +  ( KWIGR*BETAZ0*OMEGA( Jl ) *AMAG( Jl ) ) **2/ 

1  ( 4 . *NU*KWIGL**2*(GAM0  -1.)) 

CONTINUE 

PLEE (1,1)  =  EWAV0/FILL 

INITIALIZE  PHASES  FOR  THE  FIRST  MODE 

DO  80  J=1 , NPART 

U0 ( J , 5 )  =  KPOD(l)  -  OMEGA ( 1 ) 

U0T(J,1)  =  KPOD(l)  -  OMEGA ( 1 ) 

PSI ( 1 , J , 5 )  «  PS00{lfJ)  +  FLOAT( NPART-  J ) *TAU*U0 ( J , 5 ) 

PSIT( 1 , J , 1 )  =  PS00<1,J)  +  FLOAT ( NPART- J ) *TAU* ( KPOD ( 1 ) —OMEGA ( 1 ) ) 

CONTINUE 

PLEE( 1,2)  =  1. 

PLEE (1,3)  =  1.  +  PLEE (1,1) 

INITIALIZE  AMPLITUDE  EVOLUTION  WITH  THREE  POINTS  FROM  RUNGE-KUTTA 
EWAV(l)  =  EWAV0 
DO  1000  N  =  2,4 
Nl  =  N  -1 

TIME ( N )  =  F  LO AT ( N 1 ) *  TAUT  +1. 

DO  85  NN=2 , 4 

DO  90  JP  =NPART-F+1 , NPART 
J  =  JP  +  Nl *  F 

U0T ( JP , NN )  =  KPOD ( 1 )  -  OMEGA ( 1 ) 

PSIT( 1 ,JP, NN ) =  PS00(1,J)  + FLOAT ( NPART- J  P ) *TAU*U0T ( JP , NN ) 
CONTINUE 

DO  1100  Jl=l , NMODE 
CALL  EVOLVR( Jl , Kll , 1 ) 

AT (  J 1 , 2 )  =  AT ( J 1 , 1 }  +  TAUT* K1 1 ( Jl ) /2 . 

APT( Jl , 2 )  =  Kll(Jl) 

APT( Jl , 1 )  =  Kll(Jl) 

PHIl ( Jl , 7-N )  -  PHIlT( Jl , 1 ) 

PHI  2 ( Jl , 7-N )  =  PHI 2T( Jl , 1 ) 

I F ( N  .EQ.  2)  THEN 
Fll(Jl)  =  Kll(Jl) 

APP( Jl , 5 )  =  Kll(Jl) 

END  IF 

I F ( N  .EQ.  3)  F12 { Jl )  =  Kll(Jl) 

I F { N  .EQ.  4)  Fl 3 ( Jl )  =  Kll(Jl) 

CONTINUE 

DO  1200  JP=1 , NPART-F 
CALL  F4R( JP, K41 , 1 ) 

K31 ( JP )  =  U0T(JP,1) 

PSIT(  1 , JP , 2 )  =  PSIT( 1 , JP , 1 )  +  TAUT  *  K  3 1 ( J  P ) /2 . 

U0T(  JP , 2 )  =  U0T(JP,1)  +  TAUT*K41 ( JP )/2 . 

DO  1300  Jl=l, NMODE 


I  CALL  EVOLVR (  J 1 , Kl 2 , 2 ) 

AT (  J 1 , 3  )  -  AT ( J 1 , 1 )  +  TAUT*Kl2( Jl)/2. 

1300  APT( Jl ,  3 )  -  Kl 2 ( Jl ) 

DO  1400  JP=1 , NPART-F 

I  CALL  F4R ( JP , K4 2 , 2 ) 

K32 ( JP )  -  U0T(JP,1)  +  TAUT*K41( JP)/2. 

PSIT(  1 , JP , 3 )  =  PSIT(  1 , JP , 1 )  +  TAUT*K32( JP)/2. 

11400  l*0T(  JP ,  3  )  =  UOT (  JP ,  1 )  +  TAUT*K4  2  (  JP  )/2  . 

DO  1500  Jl  =  l , NMODE 
CALL  EVOLVR (J1,K13, 3) 

AT( Jl , 4 )  =  AT ( J 1 , 1 )  +  TAUT*K13 ( Jl ) 

11500  APT ( Jl , 4 )  -  Kl 3 ( Jl ) 

j  DO  1600  JP=1 , NPART-F 

CALL  F4R ( JP , K4  3 , 3 ) 

K33 ( JP )  -  U0T( JP , 1 )  +  TAUT*K42 ( JP )/2  . 

PSIT( 1 , JP , 4 )  =  PSIT( 1 , JP , 1 )  +  TAUT*K3  3 ( JP ) 

11600  U0T( JP , 4 )  =  U0T( JP, 1 )  +  TAUT  *  K  4  3 ( J  P ) 

EWAV(N)  =  0. 

)  DO  1900  Jl=l , NMODE 

I  CALL  EVOLVR ( J 1 , K 1 4 , 4 ) 

A ( Jl , 6-N ) =AT ( Jl , 1 ) +TAUT* (Kll(Jl)+2.*Kl2(Jl)+2.*K13(Jl)+ 

1  K14 ( Jl ) )/6 . 

AMAG(Jl)  =  CABS (  A(J1,6-N)  ) 

PLOTA( N  r  Jl  )  =  AMAG(Jl) 

APP ( Jl , 6-N )  =  K14(J1) 

APT ( Jl , 1 )  =  APP ( Jl , 6-N ) 

AT ( J 1 , 1 )  =  A( Jl,6-N) 

FREQ( N , Jl )  =  AIMAG(  APP ( Jl , 6-N ) /A( Jl , 6-N )  ) 

PLAR ( N , Jl )  =  REAL (  A(Jl,6-N)  ) 

PLAI ( N , Jl )  =  AIMAG(  A(Jl,6-N)  ) 

EWAV(N)  =  EWAV(N)  +  KWIGR**2* ( ( BETAZ 0* OMEGA ( Jl ) *AMAG( Jl ) )  **2  + 
1  KPOD( Jl)**2*(PHIlT(Jl,l)**2+PHl2T(Jl,l)**2)/4. )/(4.*NU* 

1  (GAM0-1 . ) *KWIGL**2 ) 

1900  GROWTH ( N , J 1 )  *  REAL (  APP ( Jl , 6-N )/A ( Jl , 6-N )  ) 

GROE(N)  =  (EWAV(N)  -EWAV0 )/( TAUT*EWAV0 ) 

EWAVO  -  EWAV(N) 

PLEE ( N , 1 )  =  EWAV ( N ) /F ILL 
PLEE ( N , 2 )  =  0. 

DO  1950  JP  =1, NPART-F 
CALL  F4R( JP,K44,4) 

K34 ( JP )  =  U0T( JP , 1 )  +  TAUT*K4  3 { JP ) 

U0 ( JP , 6-N )  =  U0T( JP , 1 ) +TAUT* (K41(JP)+2.*K42(JP)+2.*K43(JP) 

1  +K44 ( JP )  )/6 . 

UOT ( JP , 1 )  =  U0 ( JP , 6-N ) 

BETAZ  =  BETAZ 0  * ( U0 ( JP , 1 )  +  OMEGA( 1 ) )/KPOD( 1 ) 

GAM  «  1 ./SQRT( 1 . -BETAZ*BETAZ-BETAW*BETAW ) 

PLEE ( N , 2 )  -  PLEE( N, 2 )  +  (GAM  -  1.) 

PSI(1, JP,6-N)  -  PSIT(1,JP,1)  +  TAUT* ( K31 ( JP )  +2.*K32(JP)  + 

1  2  .  *K3 3 ( JP )  +  K34 ( JP )  )/6 . 

1950  PSIT ( 1 , JP , 1 )  =  PSI(1,JP, 6-N) 

DO  1975  JP-NPART-F+1 ,NPART 
U0 ( JP, 6-N )  =  U0T( JP , 1 ) 

PLEE  ( N ,  2  )  ==  PLEE  ( N ,  2  )  +  (GAM0  -  1.) 

1975  PSI ( 1 , JP , 6-N )  =  PSIT(1,JP,1) 

PLEE( N, 2 )  *  PLEE ( N , 2 ) *TRISE ( N ) /( FLOAT ( NPART ) * ( GAM0  -1.)) 

PLEE( N, 3 )  =  PLEE ( N , 2 )  +  PLEE(N,1) 

1000  CONTINUE 

C***  NOW  EVOLVE  PARTICLES  AND  FIELDS  WITH  ADAMS-BASHFORTH  PREDICTOR 

C***  CORRECTOR  METHOD  USING  THE  RESULTS  OF  THE  FOUR  PREVIOUS  TIMES 
C***  AS  INITIAL  CONDITIONS 

DO  7000  N=5 , NTIMES 
Nl  =  N  -1 

TIME ( N )  =  1.  +  FLOAT ( Nl ) *TAUT 
DO  5100  Jl=l , NMODE 
5100  CALL  EVOLV( Jl,Fl4,2) 

DO  5200  JP*1 ,NPART-4*F 


1 

1 

5200 


5300 

5400 

1 

5500 

5700 


5600 

1 

1 


1 

1 

5750 


CALL  F4 ( JP+F , FOUT , 2 ) 

CALL  F4(JP+2*F , FOUT1 , 3 ) 

CALL  F4(JP+3*F , FOUT2 , 4 ) 

CALL  F4(JP+4*F, FOUT 3 , 5 ) 

U0(JP,1)  =  U0 ( JP+F , 2 )  +  TAUT* ( 55. * FOUT  -59.*FOUTl  +37.*FOUT2 
-9 . *FOUT3 )/24 . 

PSI(1,JP,1)  -  PSI ( 1 , JP+F , 2 )  +  TAUT* (55.*U0(JP+F,2)  - 
59.*U0( JP+2*F,3)+37.*U0( JP+3*F , 4 ) -9 . *U0 ( JP+4*F,5) )/24. 

TEMP (  JP )  =  19 . *FOUT  -5.*FOUTl  +  FOUT2 

TEMPI  ( JP )  «  19.*U0(JP+F,2)-5.*U0{JP+2*F,3)+U0(JP+3*F,4) 

SUM  =  0 . 

DO  5300  JP=NPART-4*F+1 , NPART-F 
CALL  F4( JP+F, FOUT, 2) 

U0(JP,1)  =  U0 ( JP+F , 2 )  +  TAUT* FOUT 

BETAZ  -  BETAZ0 * ( U0 ( Jl , 1 )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 . /SQRT( 1 .  -BETAZ* BETAZ  -BETAW*BETAW ) 

SUM  -  SUM  +  GAM  -1. 

PSI(1,JP,1)  =  PSI(1,JP+F,2) +TAUT* ( UO ( JP+F , 2 )  +U0( JP, 1 ) )/2 . 

DO  5400  JP  =  NPART-F+1 , NPART 
J  -  JP  +  Nl*F 

U0(JP,1)  =  KPOD(l)  -  OMEGA ( 1 ) 

SUM  -  SUM  +  GAM0  -1. 

PSI(1,JP,1)  -  PS00(1,J)+  FLOAT( NPART-JP ) *TAU*U0 ( JP , 1 ) 

DO  5500  Jl=l , NMODE 

A( Jl , 1 )  =  A( Jl , 2 ) +TAUT* ( 55. *F14 (Jl )-59 . *F13( Jl ) +  3  7 . *F12( Jl ) 

-  9.*F11( Jl) )/24 . 

APP ( Jl , 1 )  -  Fl 4 ( Jl ) 

ATEMP(Jl)  -  19 . *F14 ( Jl )  -  5.*F13{J1)  +  Fl2(Jl) 

CONTINUE 

DO  5800  M=1 , MAX IT 
DO  5700  Jl-1, NMODE 
CALL  EVOLV ( Jl , Fl 5 , 1 ) 

CAPP(Jl)  -  APP ( Jl , 1 ) 

CORA(Jl)  -  A( Jl , 1 ) 

A( Jl , 1 )  *  A(J1,2)  +  TAUT* (9.*F15(J1)  +  ATEMP( Jl ) )/24. 

APP ( Jl , 1 )  =  F15( Jl ) 

EWAV(N)  =  0. 

PLEE ( N , 2 )  =  0. 

DO  5600  JP  =1 ,NP ART-4 *F 
CALL  F4 ( JP , FOUT , 1 ) 

U0(JP,1)  =  U0 ( JP+F , 2 )  +TAUT* ( 9 . *  FOUT  +TEMP ( JP ) ) /2 4 . 

BETAZ  =  BETAZ 0* ( U0 ( Jl , 1 )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 . /SQRT ( 1 .  -BETAZ* BETAZ  -  BETAW*BETAW ) 

PLEE ( N , 2 )  *  PLEE (N, 2)  +  GAM  -1. 

PSI(1,JP,1)  =  PSI(1, JP+F, 2 )  +TAUT* (9.*U0(JP,1) + TEMPI ( JP ) )/24 . 
PLEE ( N , 2 )  -  ( PLEE( N , 2 )  +  SUM ) *TRI SE ( N )/( FLOAT( NPART )*( GAM0-1 .) ) 
DO  5750  Jl-1, NMODE 

TEST-ABS(AIMAG(APP(Jl,l)/A(Jl,l) ) -AIMAG ( CAPP ( Jl ) /CORA( Jl ) ) ) 

I F (  CABS (  A( Jl , 1 ) -CORA( Jl ) )/CABS ( CORA ( Jl ) )  .GT.  ERROR  .OR. 

CABS  ( A(  Jl ,  1  )  —  CORA(  Jl  )  )  /CABS  (,  CORA  (  Jl )  -  A(  Jl ,  2  )  )  .GT.  ERROR  .OR. 

TEST/ABS ( A I MAG ( CAPP ( Jl )/CORA( Jl ) ) )  .GT.  ERROR2 )  THEN 
GO  TO  5799 
ELSE 

AMAG(Jl)  =  CABS ( A( Jl , 1 ) ) 

PLOTA( N, Jl )  =  AMAG(Jl) 

GROWTH ( N , J 1 )  =  REAL (  APP( Jl , 1 )/A( Jl , 1 )  ) 

PLAR ( N , Jl )  -  REAL(  A(Jl,l)  ) 

PLAI (N, Jl )  -  AIMAG(  A(Jl,l)  ) 

FREQ( N , Jl )  -  AIMAG(  APP{ Jl , 1 )/A( Jl , 1 )  ) 

EWAV(N)  =  EWAV(N)  +  KWIGR* * 2 * ( ( BETAZ0 *OMEGA ( Jl ) *AMAG ( Jl ) ) *  * 2  + 
KPOD (Jl)**2*(PHIl(Jl,l)**2+PHl2(Jl,l)**2)/4. )/(4.*NU* 

( GAM0-1 . ) *KWIGL**2  ) 

END  IF 
CONTINUE 

GROE(N)  =  (  EWAV(N)  -EWAV0 ) /( TAUT*EWAV0 ) 

EWAV0  -  EWAV(N) 


PLEE ( N , 1 )  -  EWAV ( N )/FI LL 
I  PLEE ( N , 3 )  =  PLEE( N, 2 )  +  PLEE(N,1) 

1  GO  TO  5850 

5799  I F  ( M  .EQ.  MAXIT ) WRITE ( 6 , 1 ) N, M 

5800  CONTINUE 

| 850  DO  6200  K=4 ,1,-1 

DO  5900  Jl=l ,NMODE 
A( J1,K+1)  -  A( Jl , K ) 

900  APP( Jl , K+l )  =  APP ( Jl , K ) 

DO  6100  JP=1 , NPART 
U0 ( JP , K+l )  =  U0(JP,K) 

PSI ( 1 , JP,K+1  )  =  PSI ( 1 , JP , K ) 

100  CONTINUE 

u200  CONTINUE 

DO  6300  Jl  =  l , NNODE 
Fll(Jl)  =  F12 ( Jl ) 

F12 ( Jl )  -  Fl 3 ( Jl ) 

6300  Fl  3  (  Jl )  **  Fl  4  (  Jl ) 

7000  CONTINUE 

WRI TE ( 1 )  TIME , PLOTA , GROWTH , FREQ , EWAV , GROE , PLAR , PLAI , 

1  KPOD , OMEGA , KWIGR , NU , GAM0 , BETAW , RI SE , 

1  FILL, REF, EPS, PHASE, ERROR, BETAZ 0 , ERROR2 

WRITE ( 2 ) NWIG , NPART , F , NMODE , NPLUS , MAXIT , NSEP , NTIMES 
CLOSE ( UNIT“1 ) 

CLOSE { UNIT=2 ) 

END 

SUBROUTINE  F4R ( JP , K4 , MM ) 

REAL  BETAZ  0 , BETAZ , GAM , KPOD (20), OMEGA (20) 

REAL  PHI1T( 20,4) ,PHI2T(20,4) , PSIT( 20 , 3000 , 4 ) 

REAL  U0T( 3000,4) ,K4( 3000) 

COMPLEX  AT( 20,4) ,APT(20,4) , SUMl , SUM2 , NUM 
INTEGER  JP, MM, NMODE 

COMMON/BLKl/BETAZO , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK2/PS IT , U0T,PHIlT, PHI 2T , AT , APT 
COMMON/BLK4/ALPHA2 , ALPHA3 , NMODE 
NUM  =  CMPLX ( 0 . ,1. ) 

BETAZ  =  BETAZ  0*(U0T(JP, MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 . /SQRT ( 1 .  -BETAZ * BETAZ  -  BETAW*BETAW ) 

SUMl  -  KPOD ( 1 ) * ( PHI 2T ( 1 , MM ) *COS ( PSIT ( 1 , JP , MM ) ) - 
1  PHI IT ( 1 , MM )*SIN(PSIT{ 1, JP,MM) ) ) 

SUM2  =  ( KPOD ( 1 ) — BETAZ 0  * BETAZ  * OMEGA ( 1 ) ) * ( CONJG ( AT ( 1 , MM )  ) * 

1  CEXP ( NUM * PSIT ( 1 , JP , MM ) ) +AT( 1 , MM ) *CEXP ( -NUM*PS IT ( 1 , JP , MM ) ) ) 

1  -NUM * B ETAZ 0  * BETAZ* KPOD ( 1 )*( CONJG ( APT ( 1 , MM ) )*CEXP( 

1  NUM*PSIT( 1 , JP , MM ) ) -APT( 1 , MM ) *CEXP ( -NUM*PSIT( 1 , JP , MM ) ) ) 

DO  100  Jl=2, NMODE 

SUMl  =  SUMl  +  KPOD ( Jl ) * ( PHI 2T ( Jl , MM )*COS(PSIT(Jl,JP, MM ) ) 

1  -PHIlT(Jl , MM )*SIN(PSIT(Jl,JP, MM ) ) ) 

j  SUM2  =  SUM2  + ( KPOD( Jl ) -BETAZ 0 * BETAZ * OMEGA ( Jl ) ) * ( 

I  1  CONJG ( AT ( J 1 , MM ) ) *CEXP ( NUM* PS I T( Jl , JP , MM ) ) +AT ( Jl , MM ) * 

1  CEXP ( -NUM*PSIT ( Jl , JP , MM ) ) ) -NUM * BETAZ 0* BETAZ* KPOD ( Jl ) * 

I  1  ( CONJG ( APT ( Jl , MM ) ) *CEXP ( NUM*PSIT( Jl , JP , MM ) ) -APT ( Jl , MM ) * 

1  CEXP ( — NUM*PSIT( Jl , JP , MM ) ) ) 

100  CONTINUE 

K4(JP)  «  ALPHA2 *REAL( SUMl )*(1. -BETAZ *BETAZ) /GAM  + 

I  1  ALPHA3*REAL( SUM2 )/( GAM*GAM ) 

RETURN 

END 

SUBROUTINE  F4 ( JP , FOUT , MM ) 

REAL  BETAZ, BETAZ 0, GAM, KPOD (20) ,OMEGA(20) 

REAL  PSI (20, 3000,5) ,U0(  3000,5) ,PHIl( 20 ,5) , PHI 2 ( 20 ,5) , FOUT 
COMPLEX  A( 20 , 5 ) ,APP( 20 , 5 ) , SUMl , SUM2 , NUM 
INTEGER  JP, MM, NMODE 

COMMON/BLKl/BETAZO , GAM0 , BETAW, KPOD, OMEGA 

ICOMMON/BLK5/PSI ,U0 , PHIl , PHI 2 , A, APP 
COMMON/BLK4/ALPHA2 , ALPHA3 , NMODE 


NUM  =  CMPLX ( 0  .  ,  1 .  ) 

BETAZ  -  BETAZO* ( UO ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 ,/SQRT( 1 . -BETAZ *BETAZ  -  BETAW*BETAW) 

SUMl  =KPOD(l)*(PHl2(l,MM)*COS(PSI(l,JP,MM) ) 

1  -PHI 1 ( 1 , MM )*SIN(PSI(1,JP, MM ) ) ) 

SUM2  = ( KPOD( 1 ) -BETAZO* BETAZ *OMEGA{ 1 ) ) * ( CONJG ( A( 1 , MM ) )* 

1  CEXP ( NUM* PS I ( 1 , JP , MM ) ) +A( 1 , MM ) *CEXP ( -NUM* PS I (  1 , JP , MM ) ) ) 

1  -NUM * BETAZ 0  * BETAZ  *  KPOD ( 1 )*{ CONJG ( APP ( 1 , MM ) ) *CEXP(NUM* 

1  PSI ( 1 , JP , MM )  ) -APP ( 1 , MM ) *CEXP ( -NUM*  PS I ( 1 , JP , MM ) ) ) 

DO  100  Jl=2 , NMODE 

SUMl=SUMl  +KPOD ( Jl ) *(PHI2( Jl,MM) *COS(PSI( Jl, JP,MM) ) 

1  -PHIl( Jl,MM)*SIN(PSI( Jl, JP,MM) ) ) 

SUM2**SUM2+  (  KPOD(  Jl )  -BETAZO* BETAZ *OMEGA(  J 1 )  )*( CONJG ( A( Jl , MM )  ) 

1  *CEXP ( NUM*  PS I ( Jl , JP , MM ) ) +A( Jl , MM ) *CEXP ( -NUM*PSI ( Jl , JP , MM ) ) ) 

1  -NUM*BETAZO*BETAZ*KPOD( Jl )*( CONJG ( APP ( Jl , MM ) )*CEXP(NUM* 

1  PSI( Jl, JP,MM) ) -APP ( Jl , MM ) *CEXP ( -NUM*  PSI ( Jl , JP , MM ) ) ) 

100  CONTINUE 

FOUT  =  ALPHA2*REAL( SUMl )*( 1 .-BETAZ *BETAZ ) /GAM  + 

1  ALPHA3*REAL( SUM2 )/( GAM*GAM) 

RETURN 

END 

SUBROUTINE  EVOLVR( J 1 , Kl , MM ) 

REAL  BETAZO, BETAZ, GAM, KPOD (20) .OMEGA (20) 

REAL  PHI IT ( 20, 4 ),PHI2T(20,4),PSIT(20,3000,4) ,U0T( 3000,4) 

REAL  TIME( 5000 ) 

COMPLEX  AT( 20,4) , APT( 20,4) ,Kl(20) , DUMl , NUM , ALPHAl 
INTEGER  J 1 , MM , N , F , NPART 

COMMON/BLKl /BETAZO , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK2/PSIT, U0T, PHI IT, PHI2T, AT, APT 

COMMON/BLK 3/TIME , NPART , N , RISE , TAU , BETAl , BETA2 , F , ALPHAl 
NUM  =  CMPLX ( 0 . , 1 . ) 

DUMl  =  CMPLX ( 0 . ,0. ) 

DUM2  =  0. 

DUM3  =  0. 

DO  100  JP  =1, NPART 
J  =  JP  +  ( N- 1 ) *  F 

TRISE  «  1.  -  EXP ( -FLOAT ( J— 1 ) *TAU/RI SE ) 

BETAZ  -  BETAZ 0 * ( U0T ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 ./SQRT( 1 .-BETAZ *BETAZ  -  BETAW*BETAW ) 

PSIT(Jl,JP, MM ) =KPOD ( Jl ) *( PSIT( 1 , JP , MM ) +OMEGA( 1 ) *TIME(N) )/KPOD( 1 ) 
1  -  OMEGA ( Jl ) *TIME(N) 

DUMl=  DUMl  +  CEXP ( NUM*PSIT ( Jl , JP , MM ) ) *TAU*TRISE*GAM0/GAM 
DUM2  =  DUM2  +  COS ( PSIT( Jl , JP , MM ) ) *TAU*TRISE 
100  DUM3  =  DUM3  +  SIN(PSIT(Jl,JP, MM ) ) * TAU* TRISE 

Kl(Jl)  =  ALPHAl * AT ( J 1 , MM )  +  BETAl *DUM 1/OMEGA ( Jl ) 

PHI IT ( Jl , MM )  =  -  BETA2  *DUM2/KPOD ( Jl ) *  *  2 
PHI2T( Jl , MM )  =  -  BETA2 *DUM3/KPOD( Jl ) **2 
RETURN 
END 

SUBROUTINE  EVOLV( Jl , Kl , MM ) 

REAL  BETAZO, BETAZ, GAM, KPOD( 20) ,OMEGA(20) 

REAL  PHI 1(20, 5) , PHI 2 (20, 5) , PS  1(20, 3000, 5), U0 (3000, 5) 

REAL  TIME(  5000  ) 

COMPLEX  A( 20, 5) ,APP( 20 , 5) ,Kl( 20 ) , DUMl .ALPHAl , NUM 
INTEGER  Jl , MM, N, F, NPART 

COMMON/BLKl/BETAZO , GAM0 , BETAW, KPOD , OMEGA 
COMMON/BLK 5/P SI , U0 , PHI 1 , PHI  2, A, APP 

COMMON/B  LK  3 /T I ME , NPART , N , RISE , TAU , BETAl , BETA2 , F , ALPHAl 
NUM  =  CMPLX ( 0 . ,1. ) 

DUMl  =  CMPLX ( 0 . ,0. ) 

DUM2  -  0. 

DUM3  =  0. 

DO  100  JP  -1, NPART 
J  =  JP  +  ( N-l ) *  F 

TRISE  -  1.  -  EXP( -FLOAT! J-l )*TAU/RISE) 

BETAZ  =  BETAZO* (U0(JP, MM)  +  OMEGA ( 1 )  )/KPOD(l) 


GAM  -  1 . /SQRT ( 1 . -BETAZ *BETAZ  -  BETAW*BETAW ) 

PSI( Jl, JP,MM)«KPOD( J1)*(PSI(1, JP,MM)+OMEGA( 1 ) *TIME(N) )/KPOD(l) 
-  OMEGA ( Jl ) *TIME( N ) 

DUMl  =  DUMl  +  CEXP( NUM*PSI ( Jl , JP , MM ) ) *TAU*TRISE*GAMO/GAM 
DUM2  -  DUM2  +  COS { PSI ( J 1 , JP , MM ) ) *TAU*TRI SE 
DUM3  =  DUM3  +  SIN ( PSI ( Jl , JP , MM ) ) *TAU*TRI SE 
Kl(Jl)  »  ALPHAl *A( Jl , MM )  +  BETAl *DUMl/OMEGA ( Jl ) 

PHI 1 ( Jl , MM )  -  -  BETA2 *DUM2/KPOD ( Jl ) *  * 2 
PHI 2 ( Jl , MM )  =  -  BETA2 *DUM3/KPOD{ Jl ) **2 
RETURN 
END 
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lie  __VC$DRBl : [ MARABLE } CRAYWIG . CFT ; 1  (4036,4,0),  last  revised  on  21-MAR-1988  14: 
i,  is  a  37  block  sequential  file  owned  by  UIC  [MARABLE].  The  records  are  varia 
ble  length  with  implied  (CR)  carriage  control.  The  longest  record  is  72  bytes. 

.  )b  CRAYWIG  (1997)  queued  to  LN03QUE  on  21-MAR-1988  14:09  by  user  MARABLE,  UIC 
l MARABLE ] ,  under  account  4790  at  priority  100,  started  on  printer  LTA7 :  on 
21-MAR-1988  14:09  from  queue  VC_LN03A. 
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C***  CODE  TO  EVALUATE  TEMPORAL  EVOLUTION  OF  THE  SPRECTRA 

C*  *  *  OF  UNSTABLE  MG^ES  IN  A  HELICAL  WIGGLER  FREE  ELECTRON  LASER 

C***  DELETION  OF  FIRST  TRANSIT  TIME 

C***  FIELD  AND  PARTICLE  EQUATIONS  ARE  EVOLVED  BY  ADAMS-BASHFORTH 

C***  METHOD  WITH  INI TAL I Z AT I ON  BY  RUNGE-KUTTA  METHOD 

C***  REFORMULATION  OF  THE  PARTICLE  PHASE  3/13 
C***  CONVERSION  TO  CRAY  FORTRAN 

REAL  BETAZO , BETAO , KPOD( 20 ) ,OMEGA(20) , BETAZ , GAM 
REAL  CTHET( 20 ) , PSI ( 20, 3000, 5) ,U0( 3000 , 5) 

REAL  TEMP (3000) ,TEMPl(3000) ,ATEMP(20) ,ATEMPl(20) ,THETA(20,5) 

REAL  TP (20,5) , TIME ( 5000 ) ,PLOTA( 5000,20) , GROWTH ( 5000 , 20 ) 

REAL  FREQ (5000,20) ,EWAV(5000) ,GROE(5000) ,A(20,5) ,APP(20,5) 

REAL  PHI 1(20, 5) , PHI 2 (20, 5) , PHI1T( 20 , 4 ) , PHI2T( 20 , 4 ) ,KWIGL 
1  , CORA( 20 ) , PLAI ( 5000,20) 

REAL  KWIGR , NU , NUR , NUI , FI LL , PLEE ( 5 000 , 3 ) ,PLAR( 5000,20) 

REAL  PSIT(20,3000,4) ,U0T( 3000,4) , THETAT (20,4) ,TPT(20,4) ,AT(20,4) 
REAL  APT( 20,4) , Kl 1(20), Kl2(20),Kl 3(20), Kl 4(20), K2 1(20), K22  (  2  0 ) 
REAL  K23(20) ,K24(20) ,K31( 3000) ,K32( 3000) ,K33( 3000) 

REAL  K34( 3000) ,K41( 3000) , K42 ( 3000 ) , K4 3 ( 3000 ) , K4 4 ( 3000 ) 

REAL  Fl 1 ( 20 ) ,F12(20) , Fl 3 ( 20 ) , Fl 4 ( 20  )  , Fl 5 ( 20 ) , F21 ( 20 ) ,F22(20) 

REAL  F23(20) ,F24(20) ,F25(20) 

INTEGFR  F , MM , J  P , J , J 1 , K , MAX I T 

COMMON/BLKl /BETAZO , GAM0 , BETAW , KPOD , OMEGA 

COMMON/BLK2/PS IT , U0T,PHIlT,PHl2T, AT , THETAT , APT , TPT 

COMMON/BLK 3/TIME , NPART , N , RI SE , TAU , BETAl , BETA2 , F , NUI , NUR 

COMMON/BLK4/  ALPHAl , ALPHA2 , NMODE 

COMMON/BLK5/PSI ,U0 , PHIl , PHI  2 , A, THETA , APP , TP 

PARAMETER  ( PI=3 . 1415926535 ) 

PS00(J1,J)  =  -OMEGA( Jl ) *FLOAT( J-l ) *TAU 

TRISE(N)  =  1.  -RISE* ( EXP ( (-TIME(N)+1. ) /RISE )  -1.) 

OPEN ( UNIT= 1 , FI LE= ' FT01 ' , STATUS= ' NEW ' ) 

OPEN ( UNIT=2 ,FILE='FT02'  , STATUS= ' NEW '  ) 

1  FORMAT (  '  THE  PRED.-CORR.  METHD  FAILED  TO  CONVERGE  ON  STEP',2X, 

1  14,'  AFTER  '  , I  4 , 2X , '  INTERATIONS'  ) 

WRITE ( 6 , 10 ) 

10  FORMAT (  'INPUT  NO.  OF  PART., NO.  OF  I TERAT . , GAM , BETaWIG , KWIGR 

1  , BUDKER,NWIG, EPS, PHASE, RISE,NPLUS, NMODE, NSEP, REF, F, FILL 

1  , MAXIT, ERROR, ERROR2 ' ) 

READ ( 5, * ) NPART, NTIMES,GAM0, BETAW, KWIGR, NU,NWIG, EPS, PHASE 
1  , RISE , NPLUS , NMODE , NSEP , REF , F , FI LL , MAXIT , ERROR , ERROR2 

20  FORMAT ( '  INPUT  DATA : NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , NWIG , 

1  EPS , PHASE , RI SE , NPLUS , NMODE , NSEP , REF , F , FI LL , MAXIT , ERROR , ERROR2 ' ) 
WRITE ( 6 , 20 ) 

WRITE ( 6 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , NWIG , EPS , PHASE , 

1  R I SE, NPLUS , NMODE, NSEP, REF, F, FILL, MAXIT, ERROR, ERROR2 
KWIGL  =  2 . *  FLOAT ( NWIG ) *P I 
BETAO  =  SQRT(1.0  -  1 . 0/( GAM0*GAM0 ) ) 

BETAZO  =  SQRT( BETA0*BETA0  -  BETAW*BETAW ) 

NOPT  =  NINT( 2 . *FLOAT( NWIG ) *BETAZ0/( 1 . -BETAZO ) )  +NPLUS 
OMEGA ( 1 )  =  ( FLOAT ( NOPT ) *PI ) /BETAZ  0 
KPOD ( 1 )  =  KWIGL  +  BETAZ0*OMEGA( 1 ) 

DO  50  J=2, (NMODE-1 )/2  +  1 

OMEGA ( J )  =  FLOAT( NOPT+ (J-l ) *NSEP ) *PI/BETAZ0 
KPOD(J)  =  KWIGL  +BETAZ0*OMEGA( J ) 

OMEGA ( NMODE+2- J )  =FLOAT( NOPT- ( J-l ) *NSEP ) *PI/BETAZ0 
KPOD ( NMODE+2- J )  =  KWIGL  +  BETAZ 0* OMEGA ( NMODE+2- J ) 

50  CONTINUE 

BETAl  =  2 . *FILL*NU*KWIGL**2*BETA0*BETAW/( KWIGR**2*BETAZ0**3 ) 

NUR  =  ( 1 . -REF ) /BETAZ  0 

NUI  =  -4.*FILL*NU*KWIGL**2*(1. -BETAW* *2/2 . )/( GAM 0* BETAZO* 

1  KWTGR**2*BETAZO*OMEGA(  1 ) ) 

BETA2  =  8 . *NU*BETA0 *KWIGL* * 2/( BETAZO* KWIGR* *2 ) 

BETA2  =  0. 

ALPHAl  =  KPOD( 1 ) /BETAZO 

ALPHA2  =  . 5*BETAW*GAM0*KPOD( 1 )/BETAZ0**2 
C***  INITIALIZE  PHASE  AND  AMPLITUDE 


TIME(l)  =  1. 

IGROE(l)  =  0. 

TAU  =  1 . /FLOAT ( NPART  -1) 

TAUT  =  F  LOAT ( F ) *  TAU 

f***  INITIALIZE  EACH  MODE  UNDER  CONSIDERATION  AND  FIND 
!***  THE  WAVE  ENERGY  DENSITY  IN  THE  INITIAL  SPECTURM 
1  EWAVO  =  0. 

DO  60  Jl=l , NMODE 
I  PLOTA ( 1 , Jl )  =  EPS 

|  A( Jl , 5 )  =  EPS 

AT ( J 1 , 1  )  =  EPS 
THETA ( Jl , 5 )  =  PHASE 
THETAT( Jl ,  1  )  =  PHASE 

1  EWAVO  -  EWAVO  +  ( KWIGR*BETAZ0 *OMEGA( Jl ) *A( Jl , 5 ) ) **2/ 

1  (  4 . *NU*KWIGL**2*(GAM0  -1-)) 

' 0  CONTINUE 

PLEE (1,1)  =  EWAVO/FILL 

C***  INITIALIZE  PHASES  FOR  THE  FIRST  MODE 
DO  80  J=l, NPART 
U0 ( J , 5 )  =  KPOD(l)  -  OMEGA ( 1 ) 

UOT ( J , 1 )  =  KPOD(l)  -  OMEGA ( 1 ) 

PSI ( 1 , J , 5 )  =  PS00(1,J)  +  FLOAT { NPART-  J ) *TAU*U0 ( J , 5 ) 

PSIT( 1 ,  J  ,  1 )  =  PS00(1,J)  +  FLOAT ( NPART- J ) *TAU* ( KPOD ( 1 ) -OMEGA( 1 ) ) 
0  CONTINUE 

PLEE (1,2)  =  1. 

PLEE (1,3)  =  1.  +  PLEE (1,1) 

***  INITIALIZE  AMPLITUDE  EVOLUTION  WITH  THREE  POINTS  FROM  RUNGE-KUTTA 
EWAV(l)  =  EWAVO 
DO  1000  N  =  2,4 
Nl  =  N  -1 

TIME ( N )  =  FLOAT ( Nl ) *TAUT  +1. 

DO  85  NN=2 , 4 

DO  90  JP  =NPART-F+1 , NPART 
J  =  JP  +  Nl *  F 

UOT(JP.NN)  =  KPOD(l)  -  OMEGA ( 1 ) 

90  PSI T ( 1 , JP , NN ) =  PS00(1,J)  + FLOAT ( NPART- JP ) *TAU*U0T ( JP , NN ) 

5  CONTINUE 

DO  1100  Jl=l, NMODE 

CALL  EVOLVR( Jl , K21 , Kll , 1 ) 

TPT ( Jl , 2 )  =  K21 ( Jl  ) 

TPT ( Jl , 1  )  =  K21 ( Jl  ) 

THETAT( Jl , 2 )  =  THETAT( Jl , 1 )  +  TAUT*K21 ( Jl )/2 . 

AT ( J 1 , 2 )  =  AT ( J 1 , 1 )  +  TAUT* Kll ( Jl )/2 . 

I  APT ( Jl , 2 )  =  Kll(Jl) 

j  APT ( Jl , 1 )  =  Kll(Jl) 

PHIl (  Jl , 7-N )  =  PHI IT ( J 1 , 1 ) 

PHI2 ( Jl , 7-N )  =  PHI2T( Jl , 1 ) 

|  IF(N  .EQ.  2)  THEN 

'  Fll(Jl)  =  Kll(Jl) 

F21(J1)  =  K2 1 ( Jl ) 

j  TP( Jl , 5 )  =  K21 ( Jl ) 

(  APP( Jl , 5 )  =  Kll(Jl) 

FREQ ( 1 , Jl  )  =  TPT ( Jl , 1  ) 

GROWTH ( 1 , Jl )  =  APP(Jl,5)/A(Jl,5) 

END  IF 

I F ( N  .EQ.  3)  THEN 
Fl2 ( Jl )  =  Kll(Jl) 

F22 ( Jl )  =  K21 ( Jl ) 

END  IF 

I F ( N  .EQ.  4)  THEN 
Fl 3 ( Jl  )  =  Kll(Jl) 

F23 ( Jl )  -  K21 ( Jl  ) 

END  IF 

1100  CONTINUE 

DO  1200  JP=1 , NPART-F 
CALL  F4R( JP, K41 , 1  ) 


i 


IF 


K31 ( JP )  =  UOT ( JP , 1 ) 

PS  IT ( 1 , JP , 2 )  =  PSIT{ 1 , JP , 1 )  +  TAUT*K31 ( JP )/2 . 

1200  l)OT(  JP  ,  2  )  =  UOT  (  JP  ,  1  )  +  TAUT  *K41(JP)/2. 

DO  1300  Jl  =  l , NMODE 
CALL  EVOLVR (Jl,K22,Kl2,2) 

TPT ( Jl , 3 )  =  K22 ( Jl ) 

THETAT( Jl ,  3  )  =  THETAT ( Jl , 1 )  +  TAUT*K22 ( Jl )/2 . 

AT ( Jl , 3 )  =  AT ( Jl , 1 )  +  TAUT*K12 ( Jl )/2 . 

1300  APT ( J 1 , 3 )  =  Kl 2 ( J 1 ) 

DO  1400  JP  =  1 , NPART-F 
CALL  F4R( JP, K42 , 2  ) 

K32 ( JP )  =  U0T( JP , 1 )  +  TAUT  *K41(JP)/2. 

PS IT ( 1 , JP , 3  )  »  PSIT( 1 , JP , 1 )  +  TAUT* K 3 2 ( JP ) /2 . 

1400  U0T( JP , 3 )  =  UOT ( JP , 1 )  +  TAUT  *K42(JP)/2. 

DO  1500  Jl=l, NMODE 
CALL  EVOLVR (Jl/K23,K13,3) 

TPT( Jl ,  4  )  =  K23 ( Jl ) 

THETAT( Jl , 4 )  =  THETAT(Jl,l)  +  TAUT*K23(Jl) 

AT ( J 1 , 4 )  =  AT ( Jl , 1 )  +  TAUT*Kl 3 ( Jl ) 

1500  APT ( Jl , 4  )  =  Kl 3 ( Jl ) 

DO  1600  JP=1, NPART-F 
CALL  F4R( JP, K43 , 3  ) 

K3  3 ( JP )  =  UOT ( JP , 1 )  +  TAUT*K42 ( JP )/2 . 

PSIT( 1 , JP , 4  )  =  PS I T ( 1 , JP , 1 )  +  TAUT*  K 3  3 ( JP ) 

1600  U0T( JP , 4 )  =  UOT ( JP , 1 )  +  TAUT*K4  3 ( JP ) 

DO  1700  Jl=l, NMODE 
1700  CALL  EVOLVR( Jl , K24 , K14 , 4 ) 

DO  1800  JP=1, NPART-F 
CALL  F4R( JP, K44 , 4  ) 

1800  K34 ( JP )  =  UOT ( JP , 1 )  +  TAUT*K43(JP) 

EWAV(N)  =  0. 

DO  1900  Jl  =  l, NMODE 

A( Jl ,6-N)=AT( Jl , 1 )+TAUT*( Kll( Jl )+2.*Kl2(Jl)+2.*K13(Jl)+ 

1  Kl4(Jl))/6. 

PLOTA( N, Jl )  =  A( Jl , 6-N ) 

APP ( Jl , 6-N )  =  ( A ( J 1 , 6-N )  -  AT ( J 1 , 1 ) ) /TAUT 
APT ( J 1 , 1 )  =  APP ( Jl , 6— N ) 

AT ( Jl , 1 )  =  A ( J 1 , 6-N ) 

THETA ( Jl , 6-N ) =THETAT ( Jl , 1 )+TAUT* (K21(Jl)+2.*K22(Jl)+2.*K23(Jl) 

1  +  K24( Jl)  )/6. 

TP ( J 1 , 6-N )  =  ( THETA( Jl , 6-N)  -  THETAT ( J 1 , 1 ) ) /TAUT 
TPT ( Jl , 1  )  =  TP ( J 1 , 6-N ) 

THETAT ( Jl , 1 )  =  THETA( Jl , 6-N ) 

FREQ( N, Jl )  =  TP ( Jl , 6-N ) 

PLAR ( N , Jl )  =  AT( Jl , 1 ) *SIN( THETAT (  Jl , 1 )  ) 

PLAI ( N , Jl )  =  -AT { Jl , 1 ) *COS ( THETAT ( J 1 , 1 ) ) 

EWAV(N)  =  EWAV(N)  +  KWIGR**2 * ( ( BETAZO *OMEGA( Jl ) * A( Jl , 6-N ) ) *  * 2  + 
1  KPOD( Jl ) **2*( PHI IT ( Jl, 1 ) **2+PHl2T( Jl ,l)**2)/4.)/(4. *NU* 

1  ( GAMO-1 . ) *KWIGL**2 ) 

1900  GROWTH ( N , Jl )  -  APP ( Jl , 6-N ) /A ( Jl , 6-N ) 

GROE(N)  =  (EWAV(N)  -EWAVO )/( TAUT *EWAV0 ) 

EWAVO  =  EWAV(N) 

PLEE ( N , 1 )  =  EWAV(N)/FILL 
PLEE ( N , 2 )  =  0. 

DO  1950  JP  =1, NPART-F 

UO ( JP , 6-N )  =  U0T( JP , 1 ) +TAUT* (K41(JP)+2.*K42(JP)+2.*K43(JP) 

]  +K4  4 { JP )  )/6. 

U0T(JP,1)  =  UO ( JP , 6-N ) 

BETAZ  =  BETAZO  * ( UO ( JP , 1 )  +  OMEGA ( 1 )) /KPOD  (  1 ) 

GAM  =  1 ,/SQRT( 1 . -BETAZ *BETAZ-BETAW*BETAW ) 

PLEE  (  N ,  2  )  =  PLEE (  N ,  2  )  -(•  (GAM  -  1.) 

PSI ( 1 , JP , 6-N )  =  PSIT( 1 , JP , 1 )  +  TAUT* ( K  3 1 ( JP )  +2.*K32(JP)  + 

1  2 . *K33 ( JP )  +  K34 ( JP  )  )/6. 

1950  PSI T ( 1 , JP , 1 )  =  PSI ( 1 , JP,6-N) 

DO  1975  JP=NPART-F+1 ,NPART 
UO ( JP , 6-N )  =  UOT ( JP , 1 ) 


L 


i 


PLEE ( N , 2 )  =  PLEE ( N , 2 )  +  (GAMO  -  1.) 

L975  PS  I ( 1 , JP , 6— N )  =  PSIT ( 1 , JP , 1 ) 

PLEE ( N , 2 )  -  PLEE ( N , 2 ) *TRI SE ( N ) /( FLOAT ( NPART ) * ( GAMO  -1.)) 

PLEE ( N , 3 )  «  PLEE ( N , 2 )  +  PLEE(N,1) 

1000  CONTINUE 

:***  NOW  EVOLVE  PARTICLES  AND  FIELDS  WITH  ADAMS-BASHFORTH  PREDICTOR 

C***  CORRECTOR  METHOD  USING  THE  RESULTS  OF  THE  FOUR  PREVIOUS  TIMES 
C***  AS  INITIAL  CONDITIONS 

DO  7000  N=5 , NTIMES 
Nl  =  N  -1 

TIME(N)  =  1.  +  FLOAT ( Nl ) *TAUT 
DO  5100  Jl=l , NMODE 
3100  CALL  EVOLV( Jl,F24,Fl4,2) 

DO  5200  JP=1 ,NPART-4*F 
CALL  F4 ( JP+F, FOUT, 2 ) 

CALL  F4 ( JP+2*F, FOUTl , 3 ) 

CALL  F4 ( JP+3*F , FOUT2 , 4 ) 

CALL  F4(JP+4*F, FOUT 3 , 5 ) 

U0(JP,1)  =  U0 ( JP+F, 2 )  +  TAUT* { 55 .* FOUT  -59.*FOUTl  +37.*FOUT2 
1  -9 . *FOUT3 )/24 . 

PSI(1,JP,1)  =  PSI ( 1 , JP+F , 2 )  +  TAUT*(55.*U0( JP+F, 2)  - 
1  59 . *U0( JP+2*F, 3 )+37 . *U0( JP+3*F, 4 )-9 . *U0 ( JP+4 *F , 5 ) )/2 4 . 

TEMP ( JP )  =  19 . *FOUT  -5.*FOUTl  +  FOUT2 
200  TEMPI ( JP )  -  19.*U0( JP+F , 2 ) - 5 . *U0 ( JP  +  2 * F , 3 ) +U0 ( JP+3*F,4) 

SUM  =  0. 

DO  5300  JP=NPART-4 * F+ 1 , NPART-F 
CALL  F4( JP+F, FOUT, 2) 

U0(JP,1)  =  U0 ( JP+F, 2 )  +  TAUT*  FOUT 

BETAZ  =  BETAZ0* ( U0 ( Jl , 1 )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 . /SQRT ( 1 .  —BETAZ* BETAZ  -BETAW*BETAW ) 

SUM  =  SUM  +  GAM  -1. 

5300  PSI(1,JP,1)  -  PSI(1, JP+F,2)+TAUT*(U0( JP+F, 2)  +U0 ( JP , 1 ) )/2 . 

DO  5400  JP  =  NPART-F+1 , NPART 
J  -  JP  +  Nl * F 

U0(JP,1)  =  KPOD(l)  -  OMEGA ( 1 ) 

SUM  =  SUM  +  GAMO  -1. 

400  PSI(1,JP,1)  =  PS00(1,J)+  FLOAT ( NPART-JP ) *TAU*U0 ( JP , 1 ) 

DO  5500  Jl=l, NMODE 

A( Jl , 1 )  =  A( J 1 , 2 ) +TAUT* (55.*Fl4(Jl)-59.*Fl3(Jl)+37.*Fl2(Jl) 

1  -  9 . *Fll ( Jl ) )/24 . 

THETA ( Jl , 1 )  =  THETA( Jl , 2 )  +  TAUT* ( 55 . * F2 4 ( J 1 ) - 59 . * F2 3 ( J 1 ) 

1  +37.*F22(J1)  -  9 . *F21 ( Jl ) )/24 . 

APP ( Jl , 1 )  =  (  A ( Jl , 1 )  -A( Jl , 2 )  ) /TAUT 

TP ( Jl , 1 )  =  (  THETA ( Jl , 1 )  -  THETA ( J 1 , 2 )  )/TAUT 

ATEMP(Jl)  =  19 .  *F14  ( Jl )  -  5.*Fl3(Jl)  +  Fl2(J.l) 

ATEMPl  (  Jl )  =  19 . *F24 { Jl )  -  5.*F23(Jl)  +  F22(Jl) 

S500  CONTINUE 

DO  5800  M= 1 , MAX I T 
DO  5700  Jl=l, NMODE 
CALL  EVOLV (Jl,F25,Fl5,l  ) 

CORA( Jl)  =  A( Jl , 1 ) 

AT ( Jl , 1 )  =  A(J1,2)  +  TAUT* ( 9 . *  Fl 5 ( Jl  )  +  ATEMP ( J 1  )  ) /2 4 . 

CTHETlJl)  =  THETA ( Jl , 1 ) 

THETAT ( Jl , 1 )  =  THETA ( Jl , 2 )  +  TAUT* ( 9 . * F2 5 ( Jl )  +  ATEMPl ( J 1 )) /2 4 . 
TPT ( Jl , 1 )  =  ( THETAT ( Jl , 1 )  -THETA(Jl,2)  )/TAUT 
-700  APT ( Jl , 1 )  =  (  AT ( J 1 , 1 )  -  A( Jl , 2 ) ) /TAUT 
EWAV(N)  =  0. 

PLEE ( N , 2 )  =  0. 

DO  5600  JP  =1 ,NPART-4*F 
CALL  F4 ( JP, FOUT, 1 ) 

U0(JP,1)  =  U0 ( JP+F , 2 )  +TAUT* { 9 . *  FOUT  +TEMP ( JP ) ) /2 4 . 

BETAZ  =  BETAZ 0  * ( U0 ( Jl , 1 )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  l./SQRT(l.  -BETAZ* BETAZ  -  BETAW*BETAW ) 

PLEE ( N , 2 )  -  PLEE (  N ,  2  )  +  GAM  -1. 

'300  PSI(1,JP,1)  =  PSI( 1, JP+F, 2)  +TAUT* (9.*U0(JP,1) +TEMP1 ( JP ) ) /2  4 . 

PLEE ( N , 2 )  =  ( PLEE ( N , 2 )  +  SUM ) *TRI SE ( N )/( FLOAT ( NPART )*( GAMO- 1 .) ) 


DO  5750  Jl  =  l , NMODE 

I F (  ABS (  AT( Jl , 1 )-CORA( Jl) )/ABS(CORA( Jl  ) )  . GT .  ERROR  .OR. 

1  ABS ( THETAT ( Jl , 1 ) -CTHET ( Jl ) ) /ABS ( CTHET ( J 1 ) )  .GT.  ERROR2 ) THEN 

DO  5725  J2=l, NMODE 
A(  J2 , 1  )  =  AT ( J 2 , 1 ) 

THETA ( J2 , 1 )  =  THETAT ( J  2 , 1 ) 

APP( J  2 , 1 )  =  APT ( J2 , 1 ) 

5725  TP ( J2 , 1 )  =  TPT ( J2 , 1 ) 

GO  TO  5799 
ELSE 

PLOTA( N, Jl )  =  AT( Jl , 1 ) 

GROWTH ( N,  Jl  )  =  2 . * (  AT( Jl , 1 ) -A< Jl , 2 )  )/( TAUT* ( AT( Jl , 1 ) +A ( Jl , 2 ) ) 

PLAR ( N , J 1 )  =  AT{ Jl , 1 ) *SIN( THETAT ( Jl ,  1  )  ) 

PLAI(N.Jl)  =  -AT ( J 1 , 1 ) *COS ( THETAT ( J 1 , 1 ) ) 

FREQ( N , Jl )  =  TP ( Jl , 1 ) 

EWAV(N)  =  EWAV(N)  +  KWIGR**2* ( ( BETAZ 0 * OMEGA ( Jl ) *AT( Jl , 1 ) ) **2  + 

1  KPOD  (  J 1 ) *  *  2  * ( PHI1 ( Jl , 1 ) **2+PHl2 ( Jl,l)**2)/4.)/(4. *NU* 

1  ( GAM  0 - 1 .  ) *  KW I G  L  *  *  2  ) 

END  IF 

5750  CONTINUE 

GROE(N)  =  (  EWAV(N)  -EWAV0 )/( TAUT*EWAV0 ) 

EWAVO  =  EWAV(N) 

PLEE ( N , 1 )  =  EWAV ( N ) /FI LL 

PLEE ( N , 3 )  =  PLEE ( N , 2 )  +  PLEE(N,1) 

GO  TO  5850 

5799  I F (  M  .EQ.  MAXI T ) WRITE ( 6 , 1  )  N  ,  M 

5800  CONTINUE 

5850  DO  6200  K=4,l,-1 

DO  5900  Jl=l, NMODE 
A ( J 1 r  K+ 1  )  =  A ( J 1 , K ) 

APP (  Jl  , K  +  l  )  =  APP ( J 1 , K  ) 

THETA ( Jl , K+ 1 )  =  THETA ( J 1 , K ) 

5900  TP { J 1 , K  + 1  )  =  TP ( J 1 , K ) 

DO  6100  J  P=  1  , NPART 
U0( JP,K+1 )  =  U0( JP,K) 

PSI ( 1 , JP,K+1  )  =  PS  I ( 1 , JP , K ) 

6100  CONTINUE 

6200  CONTINUE 

DO  6300  Jl  =  l, NMODE 
Fll(Jl)  =  F 1 2 ( J 1  ) 

Fl2 ( Jl )  =  F13(J1) 

F13(J1)  =  F 1 4 ( J 1  ) 

F21 ( Jl  )  =  F22 ( Jl  ) 

F22 ( Jl  )  =  F2 3 ( Jl  ) 

6300  F23( Jl )  =  F24 ( Jl ) 

7000  CONTINUE 

WRITE ( 1 , * )  TIME, PLOTA, GROWTH, FREQ, EWAV, GROE, PLAR, PLAI , 

1  KPOD, OMEGA, KWIGR , NU , GAM0 , BETAW, RISE, 

1  FILL, REF , EPS , PHASE , ERROR , BETAZ0 , ERROR2 

WRITE ( 2,* )NWIG, NPART, F, NMODE, NPLUS , MAXIT , NSEP , NTIMES 
CLOSE ( UNIT=1 ) 

CLOSE ( UNIT=2 ) 

END 

SUBROUTINE  F4R( JP , K4 , MM ) 

REAL  BETAZ0, BETAZ, GAM, KPOD( 20) ,OMEGA(20) 

REAL  AT( 20,4) ,APT(20, 4 ) ,PHIlT( 20, 4 ) ,PHI2T (20,4 ),PSIT(20,3000,4) 
REAL  U0T(  3000,4  ) , THETAT ( 20,4) ,TPT( 20 , 4 ) , K4 ( 3000  ) 

INTEGER  JP, MM, NMODE 

COMMON/BLKl/BETAZO , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK2/PSIT,UOT, PHIlT, PHI2T, AT, THETAT, APT, TPT 
COMMON/BLK4/ALPHA1 , ALPHA2 , NMODE 

BETAZ  =  BETAZ0  * ( U0T ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 . /SQRT ( 1 .  -BETAZ * BETAZ  -  BETAW*BETAW ) 

SUMl  =  KPOD (1)*(PHI2T(1, MM ) *COS ( PSIT( 1 ,JP, MM ) -THETAT ( 1 , MM ) ) - 
1  PH I IT ( 1 , MM ) *SIN( PSIT( 1 , JP , MM ) -THETAT ( 1 , MM ) ) ) 

SUM2  =  ( KPOD ( 1 ) - BETAZ 0  *  BETAZ* ( OMEGA ( 1 ) +TPT ( 1 , MM ) ) ) *AT ( 1 , MM ) * 


1  S IN ( PSIT ( 1  ,  JP , MM ) -THETAT ( 1 , MM ) ) -BETAZO *BETAZ *APT ( 1 , MM ) * 

1  COS ( PSIT ( 1 , JP, MM) -THETAT ( 1 ,  MM) ) 

DO  100  J 1 =2 , NMODE 

SUMl  =  SUMl  +  KPOD( Jl)*(PHl2T( Jl , MM ) *COS ( PSIT( Jl , JP,MM)- 
1  THETAT ( Jl , MM) ) -PHI IT ( J 1 , MM ) *S IN ( PS IT ( J 1 , JP , MM ) -THETAT ( Jl , MM ) ) ) 
SUM2  =  SUM2  + ( KPOD ( J 1 ) -BETAZO  * BETAZ  * ( OMEGA ( J 1 ) +TPT (  Jl , MM )  )  )  * 

1  AT( Jl , MM) *SIN( PSIT( Jl , JP , MM ) -THETAT ( Jl , MM) ) -BETAZO *BETAZ* 

1  APT (  Jl  , MM ) *COS (PSIT(Jl,JP, MM ) —THETAT ( Jl, MM ) ) 

00  CONTINUE 

K4  (  JP )  =  ALPHAl *SUMl * ( 1 . -BETAZ * BETAZ ) /GAM  +  ALPHA2*SUM2/( GAM* 

1  GAM) 


1 


1 

1 

1 

1 

1 

1 

00 

1 

1 


1 


00 


RETURN 

END 

SUBROUTINE  F4 ( JP , FOUT , MM ) 

REAL  BETAZ, BETAZO, GAM, KPOD( 20) ,OMEGA(20) 

REAL  PSI(20,3000,5),U0(3000,5),PHI1(20,5), PHI 2 ( 20 , 5 ) ,A( 20 , 5 ) , 
APP( 20, 5 ) , TP( 20, 5 ) , THETA ( 20 , 5 ) , FOUT 
INTEGER  JP, MM, NMODE 

COMMON/BLKl /BETAZ  0 , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK5/PSI , U0 , PHI 1 , PHI  2 , A , THETA , APP , TP 
COMMON/BLK4/ALPHA1 , ALPHA2 , NMODE 

BETAZ  =  BETAZO* ( U0 ( JP , MM )  +  OMEGA( 1 )  )/KPOD(l) 

GAM  =  1 ,/SQRT( 1 . -BETAZ*BETAZ  -  BETAW* BETAW ) 

SUMl  =KPOD( 1 ) * ( PHI  2 ( 1 , MM ) *COS ( PS  I ( 1 , JP , MM ) -THETA ( 1 , MM )  ) 

—  PHI 1  (  1 , MM )*SIN(PSI(1,JP, MM ) -THETA ( 1 , MM )  )  ) 

SUM2  = ( KPOD ( 1 ) -BETAZO *BETAZ*( OMEGA ( 1 ) +TP ( 1 , MM ) ) ) *A ( 1 , MM ) * 

SIN ( PSI ( 1 , JP, MM) -THETA ( 1 , MM ) ) -BETAZ 0 *BETAZ *APP ( 1 , MM ) * 

COS ( PS I ( 1 , JP, MM) -THETA ( 1 , MM ) ) 

DO  100  Jl=2, NMODE 

SUMl  =  SUMl  +KPOD ( J 1 ) * ( PHI  2 ( Jl , MM ) *COS( PSI ( Jl ,JP, MM ) -THETA ( J 1 , MM )  ) 
-PHI 1 ( Jl , MM )*SIN(PSI(Jl,JP, MM ) -THETA ( Jl , MM ) ) ) 

SUM2  =  SUM2+ ( KPOD ( J 1 ) -BETAZO  * BETAZ  * ( OMEGA ( J1)+TP(J1, MM ) ) ) *A ( J 1 , MM ) 
*SIN( PSI ( Jl , JP , MM ) -THETA ( Jl , MM ) ) -BETAZO *BETAZ *APP ( J 1 , MM ) * 

COS (PS  I ( Jl, JP,MM)-THETA( Jl ,MM)  ) 

CONTINUE 

FOUT  =  ALPHAl *SUMl * ( 1 . -BETAZ  * BETAZ ) /GAM  +  ALPHA2  *SUM2/( GAM* 

GAM) 

RETURN 

END 

SUBROUTINE  EVOLVR ( J 1 , K2 , Kl , MM ) 

REAL  BETAZO , BETAZ , GAM, KPOD ( 20 ) , OMEGA ( 20 ) 

REAL  AT(20,4),PHI1T(20,4),PHI2T(20,4),PSIT(20,3000,4), 

U0T (  3000 , 4  ) 

REAL  K2 ( 20 ) , Kl ( 20 ) , TIME( 5000 ) , THETAT ( 20 , 4 ) , APT( 20 , 4 ) ,TPT( 20 , 4 ) 

REAL  NUI , NUR 

INTEGER  J 1 , MM , N , F , NPART 

COMMON/BLK1/BETAZ0 , GAM0 , BETAW, KPOD , OMEGA 
COMMON/BLK2/PSIT,U0T, PHI IT, PHI 2T, AT, THETAT, APT, TPT 
COMMON/BLK 3/TIME , NPART , N , RI SE , TAU , BETAl , BETA2 , F , NUI , NUR 
DUMl  =  0. 

DUM2  =  0. 

DUM3  =  0. 

DUM4  =  0. 

DO  100  JP  =1, NPART 
J  =  JP  +  ( N-l ) *F 

TRISE  =  1.  -  EXP ( - FLOAT ( J- 1 ) *TAU/R I SE ) 

BETAZ  =  BETAZ  0  * ( U0T ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 . /SQRT ( 1 . -BETAZ  * BETAZ  -  BETAW  *  B  ETAW ) 

PSIT( Jl ,JP, MM ) =KPOD ( J 1 ) * ( PSIT( 1 , JP , MM ) +OMEGA( 1 ) *TIME ( N ) )/KPOD(  1 ) 

-  OMEGA{ Jl ) *TIME(N) 

DUMl  =  DUM1+COS ( PSIT( Jl , JP , MM ) -THETAT ( Jl , MM ) ) *TAU*TRI SE*GAM0/GAM 
DUM2  =  DUM2+SIN( PSIT ( Jl,JP, MM) -THETAT ( Jl, MM) ) *TAU*TRISE*GAM0/GAM 
DUM3  =  DUM3  +  COS ( PS  I T ( Jl , JP , MM ) -THETAT ( Jl , MM ) ) *TAU*TRI SE 
DUM4  =  DUM4  +  SIN( PS IT( Jl , JP , MM ) -THETAT( Jl , MM ) ) *TAU*TRISE 
K1(J1)  =  -NUR* AT ( Jl , MM )/2 . -  BETAl *DUM2/OMEGA( Jl ) 


K2(J1)  =  -NUI/2 .  +  BETAl *DUMl/( AT ( Jl , MM ) * OMEGA ( Jl ) ) 

PHI IT ( Jl , MM )  =  -  BETA2*DUM3/KPOD( Jl ) **2 
PHI 2T ( J 1 , MM )  =  -  BETA2*DUM4/KPOD( Jl ) **2 
RETURN 
END 

SUBROUTINE  EVOLV( Jl , K2 , Kl , MM ) 

REAL  BETAZO , BETAZ , GAM , KPOD <  20 ) ,OMEGA(20) 

REAL  A( 20 , 5 ) , PHI 1 ( 20 , 5 ) , PHI  2 ( 20 , 5 ) ,PSI(20,3000,5) ,U0( 3000,5) 

REAL  K2 ( 20 ) , Kl ( 20 ) , TIME( 5000) ,THETA(20,5) ,APP(20,5) , TP ( 20,5) 

REAL  NUI,NUR 

INTEGER  J 1 , MM , N , F , NPART 

COMMON/BLKl/BETAZO , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK5/PSI ,U0 , PHIl , PHI  2 , A , THETA , APP , TP 
COMMON/BLK 3/TIME , NPART , N , RISE , TAU , BETAl , BETA2 , F , NUI , NUR 
DUMl  =  0. 

DUM2  =  0. 

DUM3  =  0. 

DUM4  =  0. 

DO  100  JP  =1, NPART 
J  =  JP  +  ( N— 1 ) *  F 

TRISE  =  1.  -  EXP ( -FLOAT ( J- 1 ) *TAU/RI SE ) 

BETAZ  =  BETAZO  * ( U0 ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  1 ./SQRT( 1 . -BETAZ*BETAZ  -  BETAW* BETAW ) 

PS  I ( Jl , JP,MM)=KPOD( Jl ) *( PS I ( 1 , J  P , MM ) +OMEGA ( 1 ) *TIME ( N )  )/KPOD( 1 ) 

-  OMEGA ( Jl)*TIME(N) 

DUMl  =  DUMl  +  COS ( PS  I { J 1 , JP , MM ) -THETA ( J 1 , MM ) ) *TAU*TRI SE*GAM0/GAM 
DUM2  =  DUM2  +  SIN(PSI(Jl,JP, MM ) -THETA ( Jl , MM ) ) *TAU*TRI SE*GAM0/GAM 
DUM3  =  DUM3  +  COS ( PSI ( J 1 , JP , MM ) -THETA ( Jl , MM ) ) *TAU*TRI SE 
DUM4  =  DUM4  +  SIN ( PS  I ( J 1 , JP , MM ) -THETA( Jl , MM ) ) *TAU*TRI SE 
Kl(Jl)  =  -NUR*A(Jl,MM)/2.-  BETAl *DUM2/OMEGA( Jl ) 

K2(J1)  =  -NUI/2.  +  BETAl *DUMl/( A( J 1 , MM ) * OMEGA ( Jl ) ) 

PHI 1 ( Jl , MM )  =  -  BETA2 *DUM3/KPOD ( Jl ) **2 
PHI 2 ( J 1 , MM )  =  -  BETA2*DUM4/KPOD( Jl ) **2 
RETURN 
END 
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"5  0  or 


C***  CODE  TO  EVALUATE  TEMPORAL  EVOLUTION  OF  THE  SPRECTRA 
i***  OF  UNSTABLE  MODES  IN  A  HELICAL  WIGGLER  FREE  ELECTRON  LASER 

!***  DELETION  OF  FIRST  TRANSIT  TIME 

C*  *  *  FIELD  AND  PARTICLE  EQUATIONS  ARE  EVOLVED  BY  ADAMS-BASHFORTH 

/'***  METHOD  WITH  INI TALI ZATION  BY  RUNGE-KUTTA  METHOD 

***  REFORMULATION  OF  THE  PARTICLE  PHASE  3/13 
***  CONVERSION  TO  CRAY  FORTRAN 

***  INCLUSION  OF  FREQUENCY  SHIFT  ERROR  CHECK  IN  ADAMS-BASHFORTH 
***  EQUATION  SOLVER 

***  REPLACE  EXPRESSION  FOR  THE  DERIVATIVES  WITH  THE  FUNCTIONAL 
***  EVALUATION  OF  THE  DIFFERENTIAL  EQUATION 
***  PLOT  DATA  AFTER  EVERY  10  CALCULATIONS 
***  OUTPUT  DATA  IF  TIME  LIMIT  IS  APPROACHED 
v.***  MODIFICATIONS  TO  PRODUCE  RESTART  DATA  8/1 

REAL  BETAZO , BETAO , KPOD{ 20 ) ,OMEGA(20) , BETAZ , GAM 
REAL  CTHET( 20 ) , PSI ( 20 , 3000 , 5 ) ,U0(3000,5) ,TIM 

REAL  TEMP( 3000 ) , TEMPI ( 3000) ,ATEMP( 20) ,ATEMPl( 20 ) , THETA ( 20,5) 
REAL  TP (20, 5) , TIME (5000) ,TPLOT(6) 

REAL  FREQ ( 5000,20) ,EWAV( 5000) ,A( 20, 5 ) ,APP(20, 5) 

REAL  PHI1(20,5),PHI2(20,5) , KWIGL , PSI I ( 600 , 6 ) , U0I  (  600 , 6 ) 

1  ,CORA(20) ,TPC(20) , PLAI ( 5000,20) 

REAL  KWIGR , NU , NUR , NUI , FI LL , PLAR ( 5000,20) 

REAL  Kll ( 20 ) , K12 ( 20 ) , K13 ( 20 ) ,K14(20) ,K21(20) ,K22(20) 

REAL  K23(20) ,K24(20) ,K31( 3000) ,K32( 3000) ,K33(3000) 

REAL  K3 4(3000), K4 1(3000), K4 2(3000), K4 3(3 000)  ,K44(  3000) 

REAL  FI 1(20), FI 2(20), Fl 3(20), FI 4(20), Fl 5(20), F21(20),F22(20) 
REAL  F23(20) , F24 ( 20 ) , F25 ( 20 ) 

INTEGER  F , MM ,JP,J,Jl,K, MAXI T , TL I M , NCOUNT , NLAST 
COMMON/BLKl/BETAZ 0 , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK3/TIM , NPART , N , RI SE , TAU , BETA1 , BETA2 , F , NUI , NUR 
COMMON/BLK4/  ALPHAl , ALPHA2 , NMODE 
COMMON/BLK5/PS I ,U0,PHI1,PHI2,A, THETA , APP , TP 
PARAMETER  (PI-3.1415926535) 

PS00(J1,J)  =  -OMEGA ( Jl ) *  FLOAT ( J-l ) *TAU 
TRISE(N)  =  1.  -RI SE* ( EXP ( ( -TIM+1 . ) /RISE )  -1.) 

OPEN ( UNIT-1 , FILE- ' FT 01 ' , FORM- ' UNFORMATTED ' , STATUS= ' NEW ' ) 

OPEN ( UNIT— 2 , FI LE= ' FT 02 ' , FORM- ' UNFORMATTED ' , STATUS- ' NEW ' ) 

OPEN( UNIT- 3 , FILE- ' FT03 ' , FORM- ' UNFORMATTED ' , STATUS- ' NEW' ) 

1  FORMAT (  '  THE  PRED.-CORR.  METHD  FAILED  TO  CONVERGE  ON  STEP',2X, 

1  17,'  AFTER  ' , 1 4 , 2X ,  '  INTERATIONS ’ ,  2X ,  ' ON  MODE  ',13) 

WRITE (  6 , 10 ) 

.0  FORMAT (  'INPUT  NO.  OF  PART., NO.  OF  ITERAT ., GAM , BETAWIG , KWIGR 

1  , BUDKER, NWIG, EPS, PHASE, RISE,NPLUS, NMODE, NSEP, REF, F, FILL 

1  , MAXIT, ERROR, ERROR2, TLIM ' ) 

READ ( 5 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , NWIG , EPS , PHASE 
1  ,RISE,NPLUS , NMODE, NSEP, REF, F, FILL, MAXIT, ERROR, ERROR2 , TLIM 
7.0  FORMAT (  '  INPUT  DATA :  NPART ,  NTIMES  ,  GAM0  ,  BETAW ,  KWIGR ,  NU , NWIG , 

1  EPS, PHASE, RISE, NPLUS, NMODE, NSEP, REF, F, FILL, MAXIT, ERROR, ERROR2 ' ) 
WRITE  (6,20) 

WRITE (6,*) NPART, NTIMES, GAM0, BETAW, KWIGR, NU, NWIG, EPS, PHASE, 

1  RISE , NPLUS , NMODE , NSEP , REF , F , FI LL , MAXIT , ERROR , ERROR2 , TLIM 

NCOUNT  =  1 

KWIGL  =  2 . *  FLOAT ( NWIG ) *PI 

BETAO  =  SQRT(1.0  -  1 . 0/( GAM0 *GAM0 ) ) 

BETAZO  =  SQRT( BETA0*BETA0  -  BETAW* BETAW ) 

NOPT  =  NINT{ 2 . * FLOAT ( NWIG ) * BETAZ 0/( 1 . -BETAZO ) )  +NPLUS 
OMEGA ( 1 )  =  ( FLOAT (NOPT)* PI ) /BETAZO 
KPOD(l)  -  KWIGL  +  BETAZ0*OMEGA( 1 ) 

DO  50  J-2 ,  ( NMODE- 1 ) /2  +  1 

OMEGA ( J )  =  FLOAT( NOPT+ (J-1)*NSEP) *PI/BETAZ0 
KPOD(J)  =  KWIGL  +BETAZ0*OMEGA( J) 

OMEGA( NMODE+2-J )  -FLOAT ( NOPT- ( J-l ) *NSEP ) *PI/BETAZ0 
KPOD( NMODE+2-J )  =  KWIGL  +  B ETAZ0* OMEGA ( NMODE+2-J ) 

50  CONTINUE 

BETAl  -  2 . *FILL*NU*KWIGL**2*BETA0*BETAW/( KWIGR* * 2 *BETAZ0 *  *  3 ) 

NUR  =  ( 1 . -REF )/BETAZ0 


c*** 


c*** 

c*** 


60 

c*** 


80 

C*** 


1100 


1200 


NUI  -  -4.*FILL*NU*KWIGL**2*(l.-BETAW**2/2. )/( GAM0*BETAZ0* 

1  KWIGR**2*BETAZO*OMEGA( 1 ) ) 

BETA2  -  8 . *NU*BETA0*KWIGL**2/(BETAZ0*KWIGR**2 ) 

BETA2  =■  0. 

ALPHAl  =  KPOD( 1 )/BETAZ0 

ALPHA2  =  . 5*BETAW*GAMO*KPOD{ 1 )/BETAZ0**2 
INITIALIZE  PHASE  AND  AMPLITUDE 
TIME ( 1 )  *  1. 

TAU  -  1 ./FLOAT (NPART  -1) 

TAUT  =  FLOAT ( F ) *TAU 

INITIALIZE  EACH  MODE  UNDER  CONSIDERATION  AND  FIND 
THE  WAVE  ENERGY  DENSITY  IN  THE  INITIAL  SPECTURM 
EWAV0  =  0. 

DO  60  Jl*l , NMODE 

PLAR( 1 , Jl )  =  EPS*SIN( PHASE) 

PLAI ( 1 , Jl )  *  -EPS*COS( PHASE) 

A  (  J 1 , 5)  *  EPS 
THETA ( Jl , 5 )  =  PHASE 

EWAV0  =  EWAV0  +  ( KWIGR*BETAZ0*OMEGA( Jl ) *A( Jl , 5 ) ) **2/ 

1  ( 4 . *NU*KWIGL**2*(GAM0  -1.)) 

CONTINUE 

INITIALIZE  PHASES  FOR  THE  FIRST  MODE 

DO  80  J=l, NPART 

U0 ( J , 5 )  =  KPOD(l)  -  OMEGA(l) 

PSI ( 1 , J , 5 )  *  PS00(1,J)  +  FLOAT( NPART-  J ) *TAU*U0 ( J , 5 ) 

I F  (  MOD ( J , 5 )  .EQ.  0)  THEN 
JJ  =  INT( J/5 ) 

U0I ( JJ , NCOUNT )  =  U0 ( J , 5  ) 

PSII(JJ, NCOUNT)  =  PSI ( 1 , J , 5 ) 

END  IF 
CONTINUE 

TPLOT( NCOUNT)  =  1. 

INITIALIZE  AMPLITUDE  EVOLUTION  WITH  THREE  POINTS  FROM  RUNGE-KUTTA 
EWAV(l)  *  EWAV0 
DO  1000  N  =  2,4 
Nl  =  N  -1 

TIM  =  FLOAT ( Nl ) *TAUT  +1. 

DO  1100  Jl-1, NMODE 

CALL  EVOLV (Jl,K21,Kll,7-N) 

TP ( Jl , 6-N )  =  K21(J1) 

TP ( Jl , 7-N )  =  K2 1 ( Jl  ) 

APP ( Jl , 6-N )  =  Kll(Jl) 

APP ( Jl , 7-N )  =  Kll(Jl) 

THETA( Jl , 6-N )  =  THETA( Jl , 7-N )  +  TAUT*K21 ( Jl )/2 . 

A( Jl , 6-N )  =  A( Jl , 7-N )  +  TAUT*Kll ( Jl )/2 . 

I F ( N  .EQ.  2)  THEN 
Fll(Jl)  =  Kll(Jl) 

F21(J1)  =  K21 ( Jl  ) 

FREQ ( 1 , Jl )  =  K2 1 ( Jl  ) 

END  IF 

I F ( N  .EQ.  3)  THEN 
F12 ( Jl )  =  Kll(Jl) 

F22 ( Jl )  =  K21 ( Jl ) 

END  IF 

I F ( N  .EQ.  4)  THEN 
Fl  3 ( Jl )  =  Kll(Jl) 

F23(Jl)  =  K21(J1) 

END  IF 
CONTINUE 

DO  1200  JP=1 , NPART- I NT ( F/2 ) 

CALL  F4R( JP+INT( F/2 ) , K41 , 7-N ) 

K31 ( JP )  =  U0( JP+INT(F/2) ,7-N) 

PSI( 1, JP,6-N)  =  PSI { 1 , JP+INT( F/2 ) , 7-N)  +  TAUT*K31 ( JP )/2 . 

U0 ( JP , 6-N )  =  U0( JP+INT(F/2) ,7-N)  +  TAUT*K4 1 ( JP ) /2 . 

DO  1250  JP=NPART- INT ( F/2 )  +1, NPART 
J  =  JP  +  (N-1)*F 


UO ( JP , 6-N )  -  KPOD(l)  -  OMEGA ( 1 ) 

.250  PSI(1,  JP,6-N)  -  PS00(1,J)  +FLOAT(  NPART-JP )  *TAU*U0  (  JP ,  6-N ) 

DO  1300  Jl-l,NMODE 
CALL  EVOLV (Jl,K22,K12, 6-N ) 
i  TP<Jl,5-N)  =  K22 ( Jl ) 

APP( Jl , 6-N)  *  K12 ( Jl ) 

THETA( Jl , 5-N )  *  THETA( Jl , 7-N )  +  TAUT*K22 ( Jl )/2 . 

1.300  A(  Jl ,  5-N )  -  A(  Jl ,  7-N )  +  TAUT*Kl2  (  Jl  )/2  . 

DO  1400  JP=1 ,NPART-INT( F/2 ) 

CALL  F4R ( JP+ INT ( F/2 ) , K4 2 , 6-N ) 

K32 ( JP )  -  U0(JP+INT< F/2) ,7-N)  +  TAUT*K41 ( JP )/2 . 

PSI ( 1 , JP , 5-N )  -  PSI ( 1 , JP+INT( F/2 ) ,7-N)  +  TAUT*K32 ( JP )/2 . 

.400  U0 ( JP , 5-N )  =  U0( JP+INT(F/2) ,7-N)  +  TAUT*K42 ( JP )/2 . 

DO  1450  JP-NPART  +1-INT( F/2 ) , NPART 
J=  JP  +  ( N-l ) *  F 

U0 ( JP , 5-N )  *  KPOD(l)  -  OMEGA ( 1 ) 

.450  PSI ( 1 , JP , 5-N )  =  PS00(1,J)  +  FLOAT( NPART-JP )*TAU*U0(JP, 6-N) 

DO  1500  Jl=l , NMODE 

CALL  EVOLV (Jl,K23,Kl3, 5-N ) 

TP( Jl ,6-N)  =  K23 ( Jl ) 

THETA( Jl , 6-N )  -  THETA( Jl , 7-N )  +  TAUT*K23(Jl) 

A( Jl , 6-N )  =  A( Jl , 7-N )  +  TAUT*Kl 3 ( Jl ) 

.500  APP ( Jl , 6-N )  =  Kl 3 ( Jl ) 

DO  1600  JP=1 ,NPART-F 
CALL  F4R(JP+F,K43,5-N) 

K33 ( JP )  =  U0 ( JP+F , 7-N )  +  TAUT*K42 ( JP )/2 . 

PSI( 1 , JP,6-N)  =  PSI(1, JP+F, 7-N)  +  TAUT*K3  3 ( JP ) 

1600  U0 ( JP , 6-N )  =  U0 ( JP+F , 7-N )  +  TAUT*K43(JP) 

DO  1650  JP=NPART+1-F, NPART 
J  =  JP  +  ( N-l  )  * F 

UO ( JP , 6-N )  =  KPOD(l)  -  OMEGA ( 1 ) 

1650  PSI ( 1 , JP , 6-N )  -  PS00(1,J)  +  FLOAT( NPART-JP ) *TAU*U0 ( JP , 6— N ) 

DO  1700  Jl=l, NMODE 
.700  CALL  EVOLV (Jl ,K24,K14, 6-N ) 

DO  1800  JP=1 , NPART— F 
CALL  F4( JP+F, FOUT, 6-N) 

K44 ( JP )  =  FOUT 

.800  K34(JP)  =  U0 ( JP+F , 7-N )  +  TAUT*K43(JP) 

DO  1900  Jl=l, NMODE 

A( Jl , 6-N)=A( Jl , 7— N ) +TAUT* { Kl 1 ( Jl )  +  2 . *Kl 2 ( Jl ) +2 . *Kl 3 ( Jl )  + 

1  K14 ( Jl ) )/6 . 

APP( Jl , 6-N )  =  K14 ( Jl ) 

THETA( Jl , 6-N ) =THETA( Jl , 7-N ) +TAUT* ( K2 1 ( Jl ) +2 . *K22( Jl )+2 . *K23( Jl ) 
I  1  +  K24( Jl ) )/6. 

I  TP ( Jl , 6-N )  =  K24 ( Jl ) 

1900  CONTINUE 

DO  1950  JP  =1 , NPART-F 

U0 ( JP , 6-N )  =  U0 ( JP , 7-N ) +TAUT* (K41(JP)+2.*K42(JP)+2.*K43(JP) 

1  +K4  4 ( JP )  )/6. 

PSI(1, JP,6-N)  =  PSI ( 1 , JP , 7-N )  +  TAUT* ( K31 ( JP )  +2.*K32(JP)  + 

1  2 . *K33 ( JP  )  +  K34 ( JP )  )/6 . 

950  CONTINUE 

DO  1975  JP=NPART-F+1 , NPART 
I  J  =  JP  +  ( N-l ) *F 

I  U0 ( JP , 6-N )  =  KPOD(l)  -  OMEGA(l) 

i 975  PSI ( 1 , JP , 6-N )  =  PS00(1,J)  +  FLOAT( NPART- JP ) *TAU*U0 ( JP , 6-N ) 

1000  CONTINUE 

:***  NOW  EVOLVE  PARTICLES  AND  FIELDS  WITH  ADAMS-BASHFORTH  PREDICTOR 

!***  CORRECTOR  METHOD  USING  THE  RESULTS  OF  THE  FOUR  PREVIOUS  TIMES 
C***  AS  INITIAL  CONDITIONS 

DO  7000  N=5 , NTIMES 
i  N1  =  N  -1 

TIM  =  1.  +  FLOAT( Nl ) *TAUT 
DO  5100  Jl=l, NMODE 

1.100  CALL  EVOLV (Jl,F24,F14,2) 

DO  5200  JP=1 ,NPART-4*F 


CALL  F4 ( JP+F , FOUT, 2 ) 

CALL  F4(JP+2*F, FOUTl , 3 ) 

CALL  F4(JP+3*F, FOUT2 , 4 ) 

CALL  F4(JP+4*F, FOUT 3 , 5 ) 

U0(JP,1)  -  UO ( JP+F, 2 )  +  TAUT* (55.* FOUT  -59.*FOUTl  +37.*FOUT2 
1  -9 . *FOUT3 )/24 . 

PSI(1,JP,1)  =  PSI ( 1 , JP+F , 2 )  +  TAUT*(55.*U0( JP+F, 2)  - 
1  59 . *U0( JP+2*F, 3 )+37 . *U0( JP+3*F, 4 )-9 . *U0( JP+4*F, 5 ) )/24 . 

TEMP(JP)  =  19 . *FOUT  -5.*FOUTl  +  FOUT2 
5200  TEMPI ( JP )  =  19 . *U0( JP+F, 2 )-5 . *U0( JP+2*F, 3 )+U0( JP+3*F, 4 ) 

DO  5300  JP=NPART-4  *  F+l , NPART-F 
CALL  F4( JP+F, FOUT, 2) 

U0(JP,1>  -  UO ( JP+F , 2 )  +  TAUT* FOUT 
5300  PSI(1,JP,1)  =  PSI(1,JP+F,2) +TAUT* ( U0 ( JP+F , 2 )  +U0 ( JP , 1 ) )/2 . 

DO  5400  JP  ««  NP ART- F+l , NPART 
J  -  JP  +  N1*F 

U0(JP,1)  ■=  KPOD(l)  -  OMEGA ( 1 ) 

5400  PSI(1,JP,1)  =  PS00(1,J)+  FLOAT( NPART- JP ) *TAU*U0 ( JP , 1 ) 

DO  5500  Jl=l , NMODE 

A( Jl , 1 )  =  A( Jl , 2 ) +TAUT* (55.*Fl4(Jl)-59.*Fl3(Jl)+37.*Fl2(Jl) 

1  -  9 . *F11 ( Jl ) )/24 . 

THETA ( Jl , 1 )  =  THETA( Jl , 2 )  +  TAUT* ( 55 . *F24 ( Jl ) -59 . *F23 ( Jl ) 

1  +37.*F22(J1)  -  9 . *F21 ( Jl ) )/24 . 

APP ( Jl , 1 )  =  Fl4 ( Jl ) 

C***  THETA  PRIME  IS  THE  AVERAGE  OF  THE  DISCRETE  AND  FUNCTIONAL 
C***  VALUES  OF  THE  DERIVATIVE 
TP ( Jl , 1 )  =  F24 ( Jl ) 

ATEMP(Jl)  =  19 . *Fl4 ( Jl )  -  5.*F13(Jl)  +  F12(J1) 

ATEMPl ( Jl )  =  19 . *F24 ( Jl )  -  5 . * F2 3 ( Jl )  +  F22(Jl) 

5500  CONTINUE 

DO  5800  M=1 , MAX IT 
DO  5700  Jl=l, NMODE 
CALL  EVOLV( Jl , F2 5 , Fl 5 , 1  ) 

CORA(Jl)  =  A( Jl , 1 ) 

A ( Jl , 1 )  =  A( Jl , 2 )  +  TAUT  *(9.*F15(J1)  +  ATEMP ( Jl )  ) /2  4 . 

CTHET(Jl)  =  THETA ( Jl , 1 ) 

THETA ( Jl , 1 )  =  THETA ( Jl , 2 )  +  TAUT* ( 9 . * F2 5 ( Jl )  +  ATEMPl ( Jl ) )/24 . 
TPC(Jl)  =  F25 ( Jl ) 

TP ( Jl , 1 )  =  F25 ( Jl ) 

5700  APP ( Jl , 1 )  =  Fl 5 ( Jl ) 

I F (  MOD ( N , 1 0 )  .EQ.  0  ) THEN 
NPLOT  =  I  NT  (  N/'l  0  )  +  1 
TIME ( NPLOT )  =  TIM 
EWAV( NPLOT)  =  0. 

END  IF 

DO  5600  JP  =1 , NPART-4 * F 
CALL  F4 ( JP , FOUT , 1 ) 

U0(JP,1)  =  U0 ( JP+F , 2 )  +TAUT* ( 9 . *  FOUT  +TEMP( JP ) )/24 . 

5600  PSI(1,JP,1)  =  PSI ( 1 , JP+F , 2 )  +TAUT* (9.*U0(JP,1) +TEMP1 ( JP ) )/24 . 

DO  5750  Jl=l , NMODE 

I F (  ABS (  A( Jl , 1 ) -CORA( Jl ) )/ABS ( CORA( Jl ) )  .GT.  ERROR  .OR. 

1  ABS ( THETA ( Jl , 1 ) -CTHET ( Jl ) ) /ABS ( CTHET ( Jl ) )  .GT.  ERROR  .OR. 

1  ABS ( TP (Jl,l)-TPC(Jl)  ) /ABS  (  TP  (  J 1 , 1 ) )  .GT.  ERROR2 ) THEN 

GO  TO  5799 
ELSE 

I F (  MOD ( N , 1 0 )  .EQ.  0 ) THEN 
NPLOT  =  I NT ( N/l 0 )  +  1 

PLAR( NPLOT , Jl )  =  A(Jl,l)*SIN( THETA ( Jl , 1 ) ) 

PLAI (NPLOT, Jl )  =  -A( Jl , 1 ) *COS ( THETA ( Jl , 1 ) ) 

FREQ( NPLOT, Jl )  =  TP(Jl,l) 

EWAV ( NPLOT )  =  EWAV(NPIOT)  +  KWIGR**2* ( ( BETAZ0*OMEGA( Jl ) * 

1  A( Jl , 1 ) ) **2  +  KPOD (Jl)**2*(PHIl(Jl,l)**2+PHl2(Jl,l)**2) 

1  /4 . )/( 4 . *NU* ( GAM0-1 . ) *KWIGL**2  ) 

END  IF 
END  IF 

5750  CONTINUE 


I  GO  TO  5850 

799  IF(M  .EQ.  MAXIT) WRITE ( 6 , 1 )N, M 

800  CONTINUE 

5850  DO  6200  K=4,l,-1 

I  DO  5900  Jl=l , NMODE 

1  A( Jl , K+l )  =  A( Jl , K ) 

APP ( Jl , K+l )  =  APP ( Jl , K ) 
j  THETA (  Jl  , K+l )  =  THETA ( Jl ,  K ) 

.900  TP  (  Jl ,  K+l )  =  TP  (  Jl ,  K  ) 

DO  6100  JP=1 ,NPART-F 
U0 ( JP+F , K+l )  =  U0(JP,K) 
i  PSI ( 1 , JP+F, K+l )  =  PSI ( 1 , JP , K ) 

|  ilOO  CONTINUE 

6200  CONTINUE 

DO  6300  Jl-1, NMODE 
!  Fll(Jl)  =  F12 ( Jl ) 

F12 ( Jl )  -  Fl 3 ( Jl ) 

Fl 3 ( Jl )  =  F14 ( Jl ) 

F21 ( Jl )  =  F22 ( Jl ) 

F22 ( Jl )  =  F23 ( Jl ) 

6300  F23 ( Jl )  =  F24 ( Jl ) 

CALL  SECOND (CPU) 

I F (  CPU  .GE.  . 97*FLOAT( TLIM)  )  GO  TO  7001 
7000  CONTINUE 

7001  NLAST  =  N 

JJ  =  0 

TPLOT( NCOUNT+1 )  =  TIM 
DO  7002  JP=5 , NPART , 5 
JJ  =  JJ  +  1 

U0I (JJ, NCOUNT+1 )  =  U0 ( JP , 2 ) 

7002  PSI I (JJ, NCOUNT+1 )  -  PSI(1,JP,2) 

WRITE ( 1 )  TIME , FREQ , EWAV , PLAR , PLAI , KPOD , OMEGA , KWIGR , 

1  NU , GAM0 , BETAW , RI SE , FILL , REF , EPS , PHASE , 

1  ERROR, BETAZ0 ,ERROR2,PSII ,U0I , TPLOT, TLIM, BETA0 

WR I TE { 2 ) NWI G , NPART , F , NMODE , NPLUS , MAXIT , NSEP , NT I MES 
1  ,NCOUNT, NLAST 

WRITE ( 3)  PSI,U0,Fll,Fl2,Fl3,F21,F22,F23 , A, APP , THETA , 

1  TP , PHI 1 , PHI  2 
CLOSE ( UNIT=1 ) 

CLOSE ( UNIT=2 ) 

CLOSE ( UNIT=3 ) 

END 

SUBROUTINE  F4R ( JP , K4 , MM ) 

REAL  BETAZ0 , BETAZ , GAM, KPOD ( 20 ) ,OMEGA(20) 

REAL  A( 20 , 5 ) , PHIl ( 20,5),PHI2(20,5),PSI(20,3000,5) 

REAL  U0( 3000, 5) , THETA ( 20 , 5 ) , TP ( 20 , 5 ) , K4 ( 3000) , APP ( 20, 5) 
INTEGER  JP, MM, NMODE 

COMMON/BLK1/BETAZ0 , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK5/PS I , U0 , PHI 1 , PH 1 2 , A , THETA , APP , TP 
COMMON/BLK4/ALPHA1 , ALPHA2 , NMODE 

BETAZ  =  BETAZ 0  * ( U0 ( JP , MM)  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  SQRT( ( 1 . + ( GAM0*BETAW ) **2 )/( 1 . -BETAZ*BETAZ )  ) 

SUMl  =  KPOD( 1 ' * ( PHI  2 ( 1 , MM ) *COS ( PS  I ( 1 , JP , MM ) -THETA ( 1 , MM ) ) - 
1  PHIl ( 1 , MM )*SIN(PSI(1,JP, MM ) -THETA ( 1 , MM ) ) ) 

SUM2  »  ( KPOD ( 1 ) -BETAZ 0  * BETAZ  * ( OMEGA ( 1 ) +TP ( 1 , MM ) ) ) *A ( 1 , MM ) * 

1  SIN( PSI( 1 , JP, MM) -THETA ( 1 , MM ) ) -BETAZ 0* BETAZ *APP ( 1 , MM ) * 

1  COS ( PSI ( 1 ; JP , MM ) -THETA( 1 , MM ) ) 

DO  100  Jl=2, NMODE 

SUMl  «  SUMl  +  KPOD( J1)*(PHI2(J1, MM ) *COS ( PSI ( Jl ,JP, MM ) - 
1  THETA( Jl ,MM) ) -PHIl ( Jl , MM ) *SIN( PSI ( Jl , JP , MM ) -THETA( Jl , MM ) ) ) 

SUM2  *  SUM2  + ( KPOD( J 1 ) -BETAZ 0* BETAZ  * ( OMEGA ( Jl ) +TP ( J 1 , MM )  )  )  * 

1  A( Jl ,MM)*SIN(PSI( Jl , JP, MM) -THETA ( Jl ,MM) ) -BETAZ 0* BETAZ* 

1  APP ( Jl , MM ) *COS (PSI(Jl,JP, MM ) -THETA ( J 1 , MM ) ) 

100  CONTINUE 

K4(JP-INT(F/2) )  =  ALPHA1*SUM1* ( 1 . -BETAZ*BETAZ )/GAM  +  ALPHA2 
1  *SUM2/(GAM*GAM) 


RETURN 

END 

SUBROUTINE  F4 ( JP , FOUT , MM ) 

REAL  BETAZ , BETAZO ,  GAM, KPOD( 20 ) ,OMEGA(20) 

REAL  PS 1(20, 3000, 5) ,U0(3000,5) , PHI 1(20, 5) ,PHI2(20,5),A(20,5), 

1  APP (20,5) , TP (20, 5) ,THETA(20,5) , FOUT 
INTEGER  JP , MM, NMODE 

COMMON/BLKl/BETAZO , GAM0 , BETAW, KPOD , OMEGA 
COMMON/BLK5/PSI , U0 , PHI 1 , PHI 2 , A , THETA , APP , TP 
COMMON/BLK4/ALPHA1 , ALPHA2 , NMODE 

BETAZ  -  BETAZO  * ( U0 ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  SQRT(  (l.+(GAM0*BETAW)**2)/(l.-BETAZ*BETAZ) ) 

SUMl  =KPOD ( 1 ) * ( PHI 2 ( 1 , MM ) *COS ( PS I ( 1 , JP , MM ) -THETA ( 1 , MM ) ) 

1  - PH 1 1 ( 1 , MM )*SIN(PSI(1,JP, MM ) - TH  E  TA ( 1 , MM ) ) ) 

SUM2  « ( KPOD ( 1 ) —BETAZO* BETAZ  * ( OMEGA ( 1 ) +TP ( 1 , MM )  )  )*A(1,MM)* 

1  SIN ( PSI ( 1 , JP, MM) -THETA ( 1 , MM ) ) -BETAZO* BETAZ  *APP ( 1 , MM ) * 

1  COS ( PSI ( 1 , JP , MM ) -THETA( 1 , MM ) ) 

DO  100  Jl=2, NMODE 

SUMl  =  SUMl  +KPOD ( J 1 ) * ( PHI  2 ( J 1 , MM )*COS(PSI( Jl , JP , MM ) -THETA ( J1 , MM )  ) 
1  -PHI 1 ( Jl , MM )*SIN(PSI(Jl,JP, MM) -THETA ( Jl , MM ) ) ) 

SUM2=SUM2+( KPOD( Jl ) -BETAZO *BETAZ* ( OMEGA ( Jl ) +TP ( Jl , MM ) ) ) *A( Jl , MM ) 
1  *SIN( PSI ( Jl , JP,MM)-THETA( Jl ,MM) ) -BETAZ 0 *BETAZ*APP ( Jl , MM ) * 

1  COS( PSI ( Jl , JP , MM ) -THETA ( Jl ,MM) ) 

100  CONTINUE 

FOUT  =  ALPHAl*SUMl* ( 1 . -BETAZ*BETAZ ) /GAM  +  ALPHA2*SUM2/( GAM* 

1  GAM) 

RETURN 

END 

SUBROUTINE  EVOLV( Jl , K2 , Kl , MM ) 

REAL  BETAZO, BETAZ, GAM, KPOD (20) ,OMEGA(20) 

REAL  A( 20, 5) , PHI 1( 20,5) , PHI 2 (20, 5) , PS 1(20, 3000, 5), U0 (3000, 5) 

REAL  K2 ( 20 ) , Kl ( 20 ) , TIM, THETA( 20 , 5 ) , APP (2 0,5) , TP (20, 5) 

REAL  NUI , NUR 

INTEGER  J 1 , MM , N , F , NPART 

COMMON/BLKl/BETAZO , GAM0 , BETAW , KPOD , OMEGA 
COMMON/BLK5/PSI ,U0 , PHI1 , PHI 2 , A, THETA , APP , TP 
COMMON/BLK3/TIM , NPART , N , RISE , TAU , BETAl , BETA2 , F,NUI , NUR 
DUMl  =  0. 

DUM2  =  0. 

DUM3  =  0. 

DUM4  =  0. 

DO  100  Jn  =1, NPART 
J  =  JP  +  ( N-l ) *F 

TRISE  =  1.  -  EXP ( -FLOAT ( J-l ) *TAU/RI SE ) 

BETAZ  =  BETAZO  * ( U0 ( JP , MM )  +  OMEGA(l)  )/KPOD(l) 

GAM  =  SQRT(  ( 1 .  +  ( GAM0* BETAW ) *  *2 )/( 1 . -BETAZ* BETAZ )  ) 

PSI( Jl, JP,MM)=KPOD( Jl)*(PSI(l, JP,MM)+OMEGA(l)*TIM)/KPOD( 1 ) 

1  -  OMEGA ( Jl ) *TIM 

DUMl  =  DUMl  +  COS (PSI(Jl,JP, MM )— THETA ( Jl , MM ) ) *TAU*TRI SE*GAM0/GAM 
DUM2  =  DUM2  +  SIN ( PSI ( Jl , JP , MM ) -THETA ( Jl , MM ) ) *TAU*TRI SE*GAM0/GAM 
DUM3  -  DUM3  +  COS ( PSI ( Jl , JP , MM ) -THETA( Jl , MM ) ) *TAU*TRISE 
100  DUM4  =  DUM4  +  SIN( PSI ( Jl , JP , MM ) -THETA( Jl , MM ) ) *TAU*TRISE 

Kl(Jl)  =  — NUR*A( Jl , MM )/2 . -  BETAl *DUM2/OMEGA( Jl ) 

K2 ( Jl )  =  -NUI/2 .  +  BETAl * DUMl/ ( A ( Jl , MM ) * OMEGA ( J 1 ) ) 

PHI 1 ( Jl , MM )  =  -  BETA2  *DUM3/KPOD ( Jl ) *  *  2 
PHI 2 ( Jl , MM )  =  -  BETA2*DUM4/KPOD( Jl ) **2 
RETURN 
END 


J.LLLLLLLLL  111111111111111111 1111111111111111111111111111 111 111 
■•LLLLLLLLL  Digital  Equipment  Corporation  -  VAX/VMS  Version  V4 . 6 
LLLLLLLLLL  1111111111111111111111111111111111111111111111111111 


LLLLLLLLLL 

LLLLLLLLLL 

LLLLLLLLLL 


I 

I 

I 

I 


M 

M 

AAA 

RRRR 

AAA 

BBBB 

L 

EEEEE 

MM 

MM 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M  M 

M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

RRRR 

A 

A 

BBBB 

L 

EEEE 

M 

M 

AAAAA 

R  R 

AAAAA 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

BBBB 

LLLLL 

EEEEE 

H 

H 

EEEEE 

L 

W 

W 

III 

GGGG 

55555 

H 

H 

E 

L 

W 

W 

I 

G 

5 

H 

H 

E 

L 

W 

W 

I 

G 

555 

HHHHH 

EEEE 

L 

W 

W 

I 

G 

5 

H 

H 

E 

L 

W  W 

w 

I 

G  GGG 

5 

H 

H 

E 

L 

WW 

WW 

I 

G  G 

5  5 

H 

H 

EEEEE 

LLLLL  W 

w 

III 

GGG 

555 

FFFFF 

000 

RRRR 

9  9 

1 

F 

0  0 

R 

R 

9  9 

11 

F 

0  0 

R 

R 

1 

FFFF  0  0 

RRRR 

9  9 

1 

F 

0  0 

R  R 

9  9 

1 

F 

0  0 

R 

R 

9 

1 

F 

OOO 

R 

R 

9 

111 

File  _VC$DRBl : [ MARABLE . TWSTAG ] HELWIG5 . FOR ; 1  (4137,1,0),  last  revised  on 
L5-FEB-1985  11:34,  is  a  17  block  sequential  file  owned  by  UIC  [MARABLE).  The 
records  are  variable  length  with  a  fixed  control  size  of  2  bytes  and  implied 
(CR)  carriage  control.  The  longest  record  is  72  bytes. 

Job  HELWIG5  (1984)  queued  to  LN03_QUE  on  21-MAR-1988  13:47  by  user  MARABLE,  UIC 
[MARABLE],  under  account  4790  at  priority  100,  started  on  printer  LTA8 :  on 
21-MAR-1988  13:47  from  queue  VC_LN03B. 
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CODE  TO  EVALUATE  TEMPORAL  EVOLUTION  OF  THE  SPRECTRA 
OF  UNSTABLE  MODES  IN  A  HELICAL  WIGGLER  FREE  ELECTRON  LASER 
DELETION  OF  FIRST  TRANSIT  TIME 

REAL  A{5000, 25), PS 1(25, 10000) ,UOLD( 10000) , TIME (5000) 

REAL  GROWTH( 5000,25) ,DUM(25) ,DUMl(25) , DUM2(25) ,DUM3(25) 

REAL  OMEGA ( 25) ,KPOD( 25) , EWAV( 5000) ,GROE( 5000) , THETAPP 
REAL  KWIGR , KWIGL , NU , NUR , NUI , APRIM ( 2 5 ) ,PHIl(25) ,PHI2(25) 

REAL  THETAP( 5000,25) , THETA (25) , FILL 
INTEGER  F 

CHARACTER* 40  XLAB , GLAB , YLAB 
CHARACTER *80  ABl , AB2 ( 25 ) , AB3 , AB4 
PARAMETER  ( PI=3 . 1415926535 ) 

PS00(J1,J)  -  -OMEGA ( Jl ) * FLOAT ( J-l ) *TAU  -  THETA ( Jl ) 

TRISE ( J ) =  1.  -  EXP( -FLOAT( J-l ) *TAU/RISE ) 

WRITE( 6,10) 

XLAB  -  '  TIME$ ' 

FORMAT {  'INPUT  NO.  OF  PART . , NO .  OF  ITERAT . , GAM , BETAWIG , KWIGR 
1  ,BUDKER,NWIG, EPS, PHASE, RISE, NPLUS , NMODE , NSEP , REF, F, FILL'  ) 

READ ( 5 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , NWIG , EPS , PHASE 
1  , RISE, NPLUS, NMODE, NSEP, REF, F, FILL 

FORMAT ( '  INPUT  DATA : NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , NWIG , 

1  EPS, PHASE, RISE, NPLUS, NMODE, NSEP, REF, F, FILL  IS:') 

WRITE ( 6 , 20 ) 

WRITE ( 6, * ) NPART, NTIMES, GAM0, BETAW, KWIGR, NU, NWIG, EPS, PHASE, 

1  RI SE, NPLUS, NMODE, NSEP, REF, F, FILL 

FORMAT ( '  WAVE  AMPLITUDE  HELWIG5  W/  R=',F5.3) 

WRITE ( AB4 , 30 )REF 

KWIGL  =  2 . * FLOAT ( NWIG ) *PI 

BETA0  =  SQRT( 1 .  -  1 . /( GAM0*GAM0 ) ) 

BETAZ0  =  SQRT(BETA0*BETA0  -  BETAW* BETAW ) 

NOPT  =  JIFIX(2.*FLOAT(NWIG)*BETAZO/(1.-BETAZO) )  +NPLUS 
OMEGA ( 1 )  =  ( FLOAT (NOPT) *PI )/BETAZ0 
KPOD(l)  =  KWIGL  +  BETAZ0*OMEGA( 1 ) 

DO  50  J  =  2 ,  ( NMODE-1 )/2  +  1 

OMEGA ( J )  =  FLOAT ( NOPT+ (J-1)*NSEP>  *PI/BETAZ0 
KPOD(J)  =  KWIGL  +BETAZ0  *OMEGA( J ) 

OMEGA( NMODE+2-J )  =FLOAT( NOPT- ( J-l ) *NSEP ) *PI/BETAZ0 
KPOD( NMODE+2— J )  =  KWIGL  +  BETAZ0*OMEGA{ NMODE+2-J ) 

CONTINUE 

WRITE ( ABl , 1 0  4 ) NPART , F , NTIMES , GAM0 , BETAW , NWIG 

FORMAT ( ' NPART” ' , I  4 , 2X , ' F= '  , 1 4 , 2X , ' NTIMES= ' ,I4,2X, ' GAM= '  ,F5.3,2X, 
1  ' BETAWIG= ' , F5 . 3 , 2X , ' NWIG= ',13) 

FORMAT (  ' N= ',14, 3X , ' NU= ' ,E10.4,3X, ' KPOD= '  , El  0.4 
1 , 3X ,  '  OMEGA=> '  ,  E10 . 4 , 2X ,  'NP='  ,14) 

WRITE(AB3, 106) KWIGL, EPS, FILL, RISE, KWIGR 

FORMAT ( ' KWIGL= ' ,F8.4,2X, 'EPS=' ,E10.4,2X, ' FILL= ' , F5.3,2X, 

1 ' RI  SE= ' ,F5.3,2X, ' KWIGR= ' , F8 . 4 ) 

BETAl  *  2 . * FI LL*NU* KWIGL* * 2 *BETA0* BETAW/ ( KWIGR**2*BETAZ0**3 ) 

NUR  «  ( 1 . -REF )/BETAZ0 

NUI  =»  -4  .  *FILL*NU*KWIGL**2*  ( 1  .-BETAW* *2/2  .  )/(  GAM0*BETAZ0* 

1  KWIGR**2*BETAZO*OMEGA( 1 ) ) 

BETA 2  =  8 . *NU*BETA0*KWIGL**2/( BETAZ0*KWIGR**2 ) 

BETA2  =  0. 

ALPHA2  =  . 5*BETAW*GAM0/BETAZ0**2 
INITIALIZE  PHASE  AND  AMPLITUDE 
TIME ( 1 )  =  1. 

GROE(l)  =  0. 

TAU  =  1 ./FLOAT (NPART  -1) 

TAUT  =  FLOAT ( F ) *TAU 
EWAV0  =  0. 

INITIALIZE  EACH  MODE  UNDER  CONSIDERATION  AND  FIND 
THE  WAVE  ENERGY  DENSITY  IN  THE  INITIAL  SPECTURM 
DO  60  Jl=l, NMODE 
APRIM(Jl)  =  0. 

THETAPP  =  0. 

THETAP ( 1 , Jl )  =  0. 
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THETA(Jl)  =PHASE 
PHIl(Jl)  =  0. 

A( 1 , Jl  )  =  EPS 
PHI2 ( Jl  )  =  0. 

EWAVO  =  EWAVO  + ( GAMO *KWIGR ) *  *2  * ( ( KPOD ( Jl ) -KWIGL ) *A( 1 , Jl )/ 

1  KWIGL ) *  *2/( 4 . *NU* ( GAMO- 1 .  ) ) 

GROWTH ( 1 , Jl  )  =  0. 

WRITE ( AB2( Jl ) , 105 ) Jl FIX (OMEGA ( Jl ) *BETAZ0/PI ) ,NU,KPOD( Jl ) 

1  , OMEGA( Jl ) , NPLUS 
60  CONTINUE 

C***  INITIALIZE  PHASES  FOR  THE  FIRST  MODE 

DO  80  J=  1 , NPART 

PSI(1,J)  =  PS00(1,J)  +  ( KPOD ( 1 ) -  OMEGA ( 1 ) ) * ( 1 . -FLOAT ( J-l ) *TAU ) 
UOLD(J)  =  KPOD ( 1 )  -  OMEGA ( 1 ) 

80  CONTINUE 

C***  EVOLVE  PHASES  OF  PARTICLES  FOR  TIMES  GREATER  THAN  ONE  TRANSIT 
C***  TIME 

EWAV { 1 )  =  EWAVO 
DO  2000  N  =  2 , NTIMES 
Nl  =  N  -1 

TIME ( N )  =  FLOAT ( Nl ) *TAUT  +1. 

C***  INITIALIZE  VARAIBLES  FOR  ENSEMBLE  AVERAGE  (SUM) 

DO  2080  Jl=l , NMODE 
DUM(Jl)  =  0. 

DUMl(Jl)  =  0. 

DUM2 ( Jl )  =  0. 

DUM3 ( Jl )  =  0. 

2080  CONTINUE 

C***  EVOLVE  PARTICLE  PHASES  FOR  PARTICLES  ALREADY  WITHIN  THE  RESONATOR 


C***  FIRST  PERFORM  BOOKEEPING  TO  LOCATE  PARTICLES  (WITH  ENTRANCE  TIMES 
C***  ASSOCIATED  WITH  PARTICLES  IN  THE  RESONATOR  ON  THE  PREVIOUS  TIME 
C***  STEP)  IN  THE  FIRST  ELEMENTS  OF  THE  PARTICLE  PHASE  ARRAY 
DO  2200  J-l,  NPART  -F 
JP  =  J  +  Nl *  F 
PSI(1,J)  =  PS I ( 1 , J+F ) 

UOLD(J)  =  UOLD( J+F ) 

DO  2210  Jl  =2, NMODE 

PSI(J1,J)  =KPOD( Jl)*(PSI(l,J)  + OMEGA ( 1 ) *TIME(N)+  THETA ( 1 ) ) / 

1  KPOD ( 1 )  -  OMEGA ( Jl)*TIME(N)  -  THETA ( Jl ) 

2210  CONTINUE 

BETAZ  =  BETAZ0  * ( UOLD ( J )  +  OMEGA ( 1 )  +THETAP ( Nl , 1 ) ) /KPOD ( 1 ) 

GAM  =  1 . /SQRT ( 1 .  -BETAZ* BETAZ  -BETAW*BETAW ) 

DUM4  =  0. 

DUM5  =  0. 

DO  2350  Jl  =  1, NMODE 

DUM4  -  DUM4  +  KPOD ( Jl ) * ( PHI 2 ( Jl ) *COS ( PSI ( Jl , J ) ) -PHI 2 ( Jl ) * 

1  SIN(PSI( Jl, J)  )  ) 

DUM5  =  DUM5  + ( KPOD ( Jl ) -BETAZ  * BETAZ 0* ( OMEGA ( J 1 ) +THETAP ( Nl , Jl ) ) ) * 

1  A(Nl,Jl)*SIN(PSI(Jl,J) ) - B ETAZ* BETAZ 0*APRIM( Jl ) *COS ( PSI ( Jl , J ) ) 
DUM(Jl)  =  DUM(Jl)  +  COS(PSI(Jl,J) )*TAU*TRISE( JP)*GAM0/GAM 
DUMl(Jl)  =  DUMl(Jl)  +  SIN(PSI( Jl, J) )*TAU*TRISE( JP)*GAM0/GAM 
DUM2 ( J 1 )  =  DUM2 ( Jl )  +  COS ( PS I ( Jl , J ) ) *TAU*TRI SE ( JP ) 

DUM3 ( Jl )  =  DUM3 ( Jl )  +  S IN ( PSI ( Jl , J ) ) *TAU*TRI SE ( JP ) 

2350  CONTINUE 

C***  EVOLVE  THE  PHASES  OF  THE  FIRST  MODE 

PSI(1,J)  =  PSI(1,J)  +  TAUT* UOLD ( J ) 

UOLD(J)  =  UOLD(J)  +  TAUT* ( ( 1 . -BETAZ* BETAZ ) *KPOD( 1 ) *DUM4/ 

1 ( GAM *B ETAZ 0 ) +ALPHA2  *KPOD ( 1 ) *DUM5/( GAM*GAM )  -THETAPP) 

C***  NOW  EVALUATE  THE  PHASES  OF  THE  OTHER  MODES  IN  TERMS  OF 

C***  THE  FIRST  MODE 

DO  2375  Jl=2, NMODE 

PSI(J1,J)  =  KPOD (Jl)*(PSI(l,J)+  OMEGA ( 1)*TIME(N)+  THETA(l))/ 

1  KPOD ( 1 )  -  OMEGA( Jl ) *TIME( N)  -  THETA ( Jl ) 

2375  CONTINUE 


2120  2200  CONTINUE 

2130  C***  NOW  COMPLETE  THE  ENSEMBLE  AVERAGE  FOR  PARTICLE  JUST  ENTERING 

2140  C***  THE  RESONATOR 

2150  DO  2600  J  =NPART  +1-F,NPART 

2160  JP  =  J  +  Nl*F 

2170  DUM4  =  0. 

2180  DUM5  =  0. 

2190  DO  2450  J 1  =  1 , NMODE 

2200  DUM4  -  DUM4  +  KPOD( Jl ) * ( PHI  2 ( Jl ) *COS ( PS00 ( Jl , JP )  ) -PHI 1 ( Jl ) * 

2210  1  SIN{PS00{ Jl, JP) ) ) 

2220  DUM5  -  DUM5  +(KPOD(Jl)-  BETAZ0*BETAZ0* ( OMEGA ( Jl ) +THETAP ( Nl , Jl ) ) ) 

22  30  1  *A(Nl, Jl ) *SIN( PS00( Jl , JP) )-BETAZ0*BETAZ0*APRIM( Jl ) *COS( PS00(  Jl 

JP)  ) 

2240  DUM(Jl)  =  DUM(Jl)  +  COS ( PS00 ( Jl , JP ) ) *TAU*TRISE(  JP ) 

2250  DUMl(Jl)  -  DUMl(Jl)  +  SIN( PS00 ( Jl , JP ) ) *TAU*TRISE ( JP ) 

2260  DUM2 ( Jl )  =  DUM2(Jl)  +  COS ( PS00 ( Jl , JP ) ) *TAU*TRISE ( JP ) 

2270  DUM3 ( Jl )  =  DUM3(Jl)  +  SIN( PS00 ( Jl , JP ) ) *TAU*TRISE ( JP ) 

2280  2450  CONTINUE 

2290  C***  EVOLVE  THE  PHASES  FOR  THE  FIRST  MODE 

2300  UOLD ( J ) =KPOD ( 1 ) -OMEGA ( 1 ) -THETAP ( Nl , 1 )  +  ( 1 . + FLOAT ( Nl ) *TAUT-FLOAT(  J 

P-1) 

2310  1  *TAU ) * ( ( 1 . -BETAZ0  *  BETAZ  0 ) *  KPOD ( 1 ) *DUM4/( GAM0  *BETAZ0 ) 

2315  1  +  ALPHA2*KPOD( 1 ) *DUM5/(GAM0*GAM0 )  -  THETAPP  ) 

2330  PSI(1,J)  =  PS00(1,JP)  + ( 1 . + FLOAT ( Nl ) * TAUT- FLOAT ( JP-1 ) *TAU ) * ( UOLD 

(  J) 

2340  1  +  KPOD ( 1 )  -  OMEGA ( 1 )  -  THETAP ( Nl , 1 ) )/2 . 

2350  C***  NOW  EVALUATE  THE  PHASES  OF  THE  OTHER  MODES  FROM  THE  FIRST  MODE 

2360  DO  2460  Jl=2, NMODE 

2  37  0  PSI ( Jl , J)  =KPOD{ Jl ) M  PS I ( 1 , J ) +OMEGA( 1 ) *TIME ( N ) + THETA ( 1 )  )/KPOD(  1 ) 

2380  1  -  OMEGA ( Jl ) *TIME(N)  -  THETA ( Jl ) 

2390  2460  CONTINUE 

2400  2600  CONTINUE 

2410  C***  NOW  EVOLVE  THE  POTENTIALS,  GROWTH  RATES  FOR  EACH  OF  THE  MODES 

2420  C***  AND  EVALUATE  THE  TOTAL  WAVE  ENERGY  DENSITY 

2430  EWAV(N)  =  0. 

2440  DO  2500  Jl=l, NMODE 

2450  A( N , Jl )  =  ( A ( Nl , Jl ) -TAUT*BETAl *DUMl ( Jl ) /OMEGA ( J 1 ) ) / ( 1 . +TAUT*NUR/ 

!.  ) 


2460  THETAP ( N , Jl )  =  -NUI/2.  +  BETAl *DUM( Jl )/( OMEGA ( Jl)*A(N,Jl) ) 

2470  PHI 1 ( Jl )  =  -  BETA2 *DUM2 ( J 1 ) /( OMEGA ( Jl ) * *2 ) 

2480  PHI2 ( Jl )  =  -  BETA2  *  DUM3 ( Jl ) /( OMEGA ( Jl ) *  *  2 ) 

2490  GROWTH ( N , J 1 ) =  2 . * ( A( N , Jl ) -A( Nl , Jl ) ) / ( TAUT* ( A ( N , Jl ) +A ( Nl , Jl ) ) ) 

2500  THETA ( Jl )  =  THETA(Jl)  +  TAUT* THE TAP ( N , Jl ) 

2510  EWAV(N)  =  EWAV(N)  +< GAM0 *KWIGR )* *2 *(( KPOD ( Jl ) -KWIGL ) *A ( N , Jl ) / 

2520  1  KWIGL ) *  *2/( 4 . *NU* ( GAM0-1 . ) ) 

2530  2500  CONTINUE 

2540  GROE(N)  =  (EWAV(N)  -  EWAV0 )/( TAUT*EWAV0 ) 

2550  EWAV0  =  EWAV(N) 

2560  THETAPP  =  2 .*( THETAP ( N , 1 ) -THETAP ( Nl , 1 ))/( THETAP ( N , 1 ) +THETAP ( Nl , 1 

)  ) 

2570  C***  REPEAT  FOR  NEXT  TIME  STEP 

2580  2000  CONTINUE 

2590  DO  3000  Jl=l, NMODE 

2600  YLAB  =  '  WAVE  AMPLITUDE$ ' 

2610  GLAB  =  '  WAVE  AMPLITUDE  VS.  TIME  MULTI$' 

2620  CALL  ANOTAT ( %REF ( XLAB ) , %  REF ( YLAB ),1,0,0,DSHL) 

2630  CALL  PWRIT( 500 , 80 , %REF ( ABl ) , 70 , 1 , 0 , 0 ) 

2640  CALL  PWRIT( 500 , 48 , %REF ( AB2 ( Jl ) ) , 72 , 1 , 0 , 0 ) 

2650  CALL  PWRIT( 500 , 16 ,  %REF ( AB 3), 75, 1,0,0) 

2660  CALL  EZXY ( TIME , A( 1 , Jl ) , NTIMES , %REF ( AB4 ) ) 

2670  YLAB  =  '  GROWTH  RATE  $' 

2680  GLAB  =  '  GROWTH  RATE  NORMALIZED  TO  TRANSIT  TIME$ ' 

2690  CALL  ANOTAT( %REF ( XLAB ), %REF ( YLAB ), 1 , 0 , 0 , DSHL ) 

2700  CALL  EZXY(TIME, GROWTH* 1 ,J1 ), NTIMES, %REF(GLAB)  ) 

2710  YLAB  =  '  FREQUENCY  SHI FT$ ' 

2720  GLAB  =  '  FREQUENCY  SHIFT  VS.  TIME  $' 
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CALL  ANOTAT ( %REF ( XLAB ) , %REF ( YLAB ) , 1 , 0 , 0 , DSHL ) 
CALL  EZXY ( TIME , THETAP { 1 , Jl ) , NTIMES , %REF ( GLAB )  ) 

3000  CONTINUE 

GLAB  =  'FIELD  ENERGY  DENSITY  VS.  TIME$ ' 

YLAB  =  '  ENERGY  DENSITY$ ' 

CALL  ANOTAT ( %REF ( XLAB ) ,%REF( YLAB) ,1 , 0, 0, DSHL) 
CALL  EZXY ( TIME, EWAV, NTIMES, %REF ( GLAB ) ) 

GLAB  =  '  RATE  OF  CHANGE  OF  ENERGY$ ' 

YLAB  =  'GROWTH  RATE$ ' 

CALL  ANOTAT ( % REF (XLAB) , %REF ( YLAB ) ,1,0,0, DSHL) 
CALL  EZXY(TIME,GROE, NTIMES, %REF ( GLAB ) ) 

STOP 

END 


pppppppppp 

pppppppppp 

pppppppppp 


“PPPPPPPPP 
ppppppppp 
. ppppppppp 


5555555555555555555555555555555555555555555555555555 
Digital  Equipment  Corporation  -  VAX/VMS  Version  V4.6 
5555555555555555555555555555555555555555555555555555 


M 

M 

AAA 

RRRR 

AAA 

BBBB 

L 

EEEEE 

MM 

MM 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M  M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

RRRR 

A 

A 

BBBB 

L 

EEEE 

M 

M 

AAAAA 

R  R 

AAAAA 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

BBBB 

LLLLL 

EEEEE 

M 

M 

AAA 

X 

X 

w 

w 

ip  (p  (Jl  1JI  fjl 

EEEEE 

M 

M 

MM 

MM 

A 

A 

X 

X 

w 

w 

T 

E 

MM 

MM 

M 

M  M 

A 

A 

X 

X 

w 

w 

T 

E 

M  M 

M 

M 

M 

A 

A 

X 

w 

w 

T 

EEEE 

M 

M 

M 

M 

AAAAA 

X 

X 

WWW 

T 

E 

M 

M 

M 

M 

A 

A 

X 

X 

ww 

ww 

T 

E 

M 

M 

M 

M 

A 

A 

X 

X 

w 

w 

T 

EEEEE 

M 

M 

FFFFF 

000 

RRRR 

7  7 

1 

F 

0 

0 

R 

R 

7  7 

11 

F 

0 

0 

R 

R 

1 

FFFF 

0 

0 

RRRR 

7  7 

1 

F 

0 

0 

R 

R 

7  7 

1 

.  # 

F 

0 

0 

R 

R 

7 

1 

•  • 

F 

000 

R 

R 

7 

111 

"ile  _VC$DRBl : [ MARABLE . FEL ] MAXWTEM . FOR ; 1  (3924,1,0),  last  revised  on 
7-OCT-1984  15:18,  is  a  34  block  sequential  file  owned  by  UIC  [MARABLE).  The 
records  are  variable  length  with  a  fixed  control  size  of  2  bytes  and  implied 
(O')  carriage  control.  The  longest  record  is  72  bytes. 

ob  MAXWTEM  (2005)  queued  to  LN03_QUE  on  21-MAR-1988  14:21  by  user  MARABLE,  UIC 
[MARABLE],  under  account  4790  at  priority  100,  started  on  printer  LTA7 :  on 
"‘l-MAR-1988  14:21  from  queue  VC_LN03A. 
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C  MAIN  PROGRAM  FOR  COMPARISON  OF  DELTA  FUNCTION  AND  STEP  FUNCTION 
C  DISPERSIONS  USING  MULLER'S  METHOD  TO  DETERMINE  THE  COMPLEX  PART 
C  OF  THE  EIGENFREQUENCIES 

INTEGER  N , NPREV, MAXIT , I , M , L , NPTS 

DOUBLE  PRECISION  EPl , EP2 , DENS , WIG , LAM , BETA, PERP , BETAl , PERPl , GAM 
DOUBLE  PRECISION  DUMl , DUM2 , DUM3 , DUM4 , X ( 9 , 50 ) , X2 ( 9 , 50 ) 

DOUBLE  PRECISION  TS , TD , LINC , LINC2 , LINC3 , Yl ( 8 ) , OMEG 
DOUBLE  PRECISION  ALPl , ALP 2 , DEL 
CHARACTER* 40  XLAB 

DOUBLE  COMPLEX  ZEROS ( 6 ) , ROOT , PREVRT( 2 0 ) 

EXTERNAL  FNl/FN2 
LOGICAL  FNREAL , FSTG 

C  FSTG  IS  A  LOGICAL  VARIABLE:  IF  TRUE  THE  FIRST  GUESS  TO  THE  ROOTS 

C  OF  THE  DISPERSION  RELATION  ARE  THE  ROOTS  FROM  THE  PREVIOUS 

C  INCREMENT  OF  K.  IF  FALSE  THE  FIRST  GUESS  IS  FROM  PREVIOUS  DATA 

C  AT  THE  SAME  K  FROM  THE  FILE  DATINI . 

COMMON  DENS , WIG, LAM, BETA, BETAl , PERP , PERPl , GAM , ALPl , ALP2 , DEL 
OPEN  (UNIT  =1,  NAME= ' OUTF I LE '  .STATUS  ='NEW') 

OPEN  (UNIT  =2,  NAME= ' DATINI '  .STATUS  ='OLD') 

FNREAL=  .FALSE. 

MAXIT  =  1000 
EPl  =  l.D-15 
EP2  =  l.D-15 
WRITE ( 6 , 1 ) 

WRITE ( 1,1) 

1  FORMAT ( '  INPUT  FOLLOWING  DATA:  DENS , WIG , LAM , LINC , LINC2 , LINC3 , BET 
1  , PERP, DEL, NPTS, FIRSTG  ') 

READ ( 5, * ) DENS, WIG, LAM, LINC, LINC2 , LINC3 .BETA, PERP , DEL , NPTS , FSTG 
READ ( 2 , * ) ( (X2(K, J) ,K=1,9) , J=1,NPTS) 

WRITE ( 1 , * )  DENS , WIG, LAM, LINC, LINC2 ,LINC3 .BETA, PERP, DEL, NPTS, FSTG 
WRITE ( 6 , * )  DENS , WIG, LAM, LINC , LINC2 , LINC3 , BETA , PERP , DEL , NPTS , FSTG 
GAMBAR  =  DSQRT ( 1 .  +  BETA*BETA  +  PERP  +WIG) 

BETA  =  BETA/GAMBAR 

PERP  =  PERP/( GAMBAR*GAMBAR ) 

WIG  =  WIG/ (GAMBAR* GAMBAR) 

DENS  =  DENS/GAMBAR 

DEL  =  DEL/GAMBAR 

BETAl  =  BETA/DSQRT( 1 .-PERP ) 

PERPl  =  PERP/ ( 1 . -PERP ) 

GAM  =  1 ./DSQRT( 1 . -PERP ) 

TD  =  PERP*GAMBAR* ( 1 .  -  . 5*WIG) 

TS  =  . 5*PERP1*GAMBAR*( 1 .  -  ,5*WIG  -  PERPl/3 . ) /GAM 
WRITE (6,3)  TD , TS 
WRITE (1,3)  TD , TS 

3  FORMAT ( '  DELTA  TEMP  =  ',Dl3.4,'  STEP  TEMP  =  '.D13.4) 

WRITE ( 6 , 2 ) 

WRITE ( 1 , 2  ) 

2  FORMAT ( '  LAM  2nd  ROOT  DLD  DTPlD  DTMlD  DLS 

1  DTP1S  DTMlS ' ) 

DO  1000  L  =  1 , NPTS 
I F (  L  .LE.  20)  LAM  =  LAM  +  LINC 

I F (  L  .GT.  20  .AND.  L  .LE.  30)  LAM  =  LAM  +  LINC2 
I F (  L  .GT.  30)  LAM  =  LAM  +  LINC3 
ALPl  =  DEL* LAM* ( 1 . -  BETA*BETA ) 

ALP2  =  1.  -  PERPl/4 .  -BETAl *  BETAl * ( 1 .  -  3.*PERPl/4.) 

X( 1 , L )  =  LAM 
N  =  6 
NPREV  =  0 


ZEROS ( 1 )  = 

DCMPLX(X2( 3 , L ) , X2 ( 2 , L ) ) 

ZEROS ( 2 )  = 

(0.  ,0. ) 

ZEROS ( 3 )  = 

(0.  ,0. ) 

ZEROS ( 4 )  = 

(0.  ,0. ) 

ZEROS ( 5 )  = 

(0.  ,0.  ) 

ZEROS ( 6 )  = 

(0. ,0. ) 

IF  (L  .GT. 

1  .AND.  FSTG)  THEN 

4900 
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5100 
5200 
5300 
5400 
5500 
5600 
5700 
5800 
5900 
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6100 
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6400 
6500 
6600 
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6800 
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7100 
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7400 
7500 
7600 
7700 
7800 
7810 
7820 
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7900 
8000 
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3400 
3500 
8600 
3700 
3800 
8900 
9000 
9100 
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9300 
9400 
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9800 
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’  9100 
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:  1500 
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’  ^620 
'630 


ZEROS ( 1 )  =  PREVRT ( 1 ) 

ZEROS ( 2 )  =  PREVRT ( 2 ) 

ZEROS ( 3 )  =  PREVRT ( 3 ) 

ZEROS ( 4 )  =  PREVRT ( 4 ) 

ZEROS ( 5 )  «  PREVRT ( 5 ) 

ZEROS ( 6 )  =  PREVRT ( 6 ) 

END  IF 

CALL  MULLER( FN1 , FNREAL , ZEROS , N , NPREV , MAXIT , EPl , EP2 , 1 ) 

C  FIND  EIGENVALUES  WITH  THE  LARGEST  GROWTH  RATES 

CALL  SORTl ( N, ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 

Yl  ( 1  )  =  DUM3 
Yl ( 2 )  =  DUM4 
X( 2 , L )  =  DUMl 
X  ( 3 , L )  =  DUM2 
PREVRT ( 1 )  =  ZEROS ( 1 ) 

PREVRT { 2 )  =  ZEROS { 2 ) 

PREVRT ( 3 )  =  ZEROS { 3 ) 

PREVRT ( 4 )  =  ZEROS ( 4 ) 

PREVRT ( 5 )  =  ZEROS { 5 ) 

PREVRT ( 6 )  =  ZEROS ( 6 ) 

C  ******  x2  IS  THE  LARGEST  GROWTH  RATE  IN  THEN  DATA  SET  FOR  FULL  DELTA 
C  ******  X3  IS  THE  REAL  FREQUENCY  CORESPONDING  TO  THE  ABOVE 
C  *******  yl  IS  THE  SECOND  LARGEST  DISTINCT  IMAGINARY  FREQUENCY 
C  *******  y2  IS  THE  REAL  FREQUENC  CORRESPONDING  TO  THE  ABOVE 
OMEG  =  DUM2  -  BETA* LAM 

Yl ( 3 )  =  LAM* *2 * ( OMEG* *2  -  DENS* ( 1 . -BETA** 2 ) ) 

Yl ( 4 )  =  DUM2**2  -{LAM  +1.)**2  -  DENS* ( 1 . -PERP/2 . ) 

Yl ( 5 )  =  DUM2**2  -(LAM-1. )**2  -  DENS*(1.  -  PERP/2.) 

NPREV  =  0 

ZEROS ( 1 )  =  DCMPLX (X2(5,L),X2(4,L)) 

ZEROS ( 2 )  =  (0. ,0.  ) 

ZEROS ( 3 )  =  ( 0  .  ,0.  ) 

ZEROS ( 4 )  =  (0.  ,0.  ) 

ZEROS ( 5  )  =  (0.,0.) 

ZEROS ( 6  )  =  (0.,0.) 

I F ( L  .GT.  1  .AND.  FSTG )  THEN 
ZEROS ( 1 )  =  PREVRT ( 7 ) 

ZEROS ( 2 )  =  PREVRT ( 8 ) 

ZEROS( 3)  =  PREVRT ( 9 ) 

ZEROS(4)  =  PREVRT (10) 

ZEROS(5)  =  PREVRT (11) 

ZEROS ( 6 )  =  PREVRT (12) 

END  IF 

CALL  MULLER ( FN2 , FNREAL, ZEROS , N , NPREV , MAXIT , EPl , EP2 , 1 ) 

CALL  SORTl ( N, ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 

X ( 4 , L )  =  DUMl 
X ( 5 , L )  =  DUM2 
PREVRT ( 7 )  =  ZEROS { 1 ) 

PREVRT ( 8 )  =  ZEROS! 2) 

PREVRT! 9)  =  ZEROS! 3) 

PREVRT (10)  =  ZEROS ( 4 ) 

PREVRT (11)  =  Z  EROS ( 5 ) 

PREVRT (12)  =  ZEROS( 6) 

C  ****  X4  IS  THE  LARGEST  GROWTH  RATE  IN  THE  DATA  SET  FOR  FULL  STEP 
C  ****  X5  IS  THE  REAL  FREQUENCY  CORESPONDING  TO  THE  ABOVE 
OMEG  =  DUM2  -  BETAl *LAM* ( 1 .  -  PERPl/4 . ) 

Yl ( 6 )  *  DENS* GAM* ( ( 1 . -BETAl **2 )-PERPl*( 1 .-3 . * BETAl* *2 )/4 .  ) 

Yl  (  6 )  =  LAM* *2* (OMEG* *2  -  Yl(6)) 

Yl  (  7 )  =  DUM2  *  *  2- ( LAM+ 1 .  )** 2-DENS* ( 1 . -PERPl/2 .  +  ( PERPl/2 .  )** 2 ) *GAM 
Yl  (  8 )  =  DUM2**2-( LAM-1 .  )** 2-DENS* ( 1 . -PERPl/2 .  +  ( PERPl/2 .  )**2)*GAM 
N  =  4 
NPREV  =  0 

ZEROS ( 1 )  =  DCMPLX ( X2(7,L) ,X2(6,L) ) 

ZEROS ( 2 )  =  (0.,0.) 

ZEROS ( 3 )  =  (0.,0.) 

ZEROS ( 4 )  =  (0.,0.) 
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12700 
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12900 
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14100 
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15200 
15300 
15400 
I*(  14  )  ) 
15500 
15600 
15700 
15800 
15900 
(14)  ) 
16000 
16100 


IF  (L  .GT.  1  .AND.  FSTG )  THEN 
ZEROS ( 1 )  =  PREVRT (13) 

ZEROS ( 2 )  =  PREVRT(14) 

ZEROS ( 3 )  =  PREVRT(15) 

ZEROS ( 4 )  =  PREVRT(16) 

END  IF 

CALL  MULLER ( FNl , FNREAL, ZEROS , N, NPREV, MAX IT , EP1 , EP2 , 2 ) 

CALL  SORTl ( N , ZEROS , DUMl ,DUM2 , DUM3 , DUM4 ) 

X( 6 , L )  =  DUMl 
X ( 7 , L )  =  BETA*LAM 
PREVRT (13)  =  ZEROS ( 1 ) 

PREVRT (14)  =  ZEROS ( 2 ) 

PREVRT (15)  =  ZEROS ( 3 ) 

PREVRT (16)  =  ZEROS ( 4 ) 

C  *****  x6  IS  THE  LARGEST  GROWTH  RATE  FOR  THE  REFERENCE  DELTA  FUNCTION 
C  *****  X7  IS  AN  APPROXIMATION  TO  THE  REAL  FREQUENCY  FOR  SMALL  DENSITIES 
NPREV  =  0 

ZEROS ( 1 )  =  DCMPLX( X2(9,L) ,X2(8,L) ) 

ZEROS ( 2 )  =  (0.,0.) 

ZEROS ( 3 )  =  (0.,0.) 

ZEROS ( 4  )  =  (0.,0.) 

IF  (L  .GT.  1  .AND.  FSTG)  THEN 
ZEROS ( 1 )  =  PREVRT( 17 ) 

ZEROS ( 2 )  =  PREVRT (18) 

ZEROS ( 3 )  =  PREVRT (19) 

ZEROS ( 4 )  =  PREVRT (20) 

END  IF 

CALL  MULLER ( FN2 , FNREAL , ZEROS , N , NPREV, MAXIT , EPl , EP2 , 2 ) 

CALL  SORTl ( N , ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 

X ( 8 , L )  =  DUMl 

X( 9 , L )  =  BETAl *  LAM 

PREVRT( 17 )  =  ZEROS ( 1 ) 

PREVRT (18)  =  ZEROS ( 2 ) 

PREVRT ( 19 )  =  ZEROS ( 3 ) 

PREVRT ( 20 )  =  ZEROS ( 4 ) 

C  ******  x8  IS  THE  LARGEST  GROWTH  RATE  FOR  THE  REFERENCE  STEP  FUNCTION 
C  *****  X9  IS  AN  APPROXIMATION  OF  THE  REAL  FREQUENCY  FOR  SMALL  DENSITIES 
WRITE (1,5)  LAM , (Yl(K) ,K=1,4) 

WRITE ( 1,5) (Yl(K) , K  =  5 , 8 ) 

WRITE (6,5)  LAM, ( Yl ( K ) , K=1 , 4 ) 

WRITE (6, 5) ( Yl ( K ) , K=5 , 8 ) 

1000  CONTINUE 

WRITE ( 1 , 6 ) 

WRITE (6,6) 

6  FORMAT ( '  LAM  FULL  DELTA  FULL  STEP  REF  DELTA 

1  REF  STEP  ' ) 

DO  2000  J  =  1 , NPTS 
WRITE (1,5) (X(K,J) , K=1 , 5 ) 

WRITE (1,5) (X(K,J)  ,K  =  6,9) 

WRITE ( 1 , * ) 

WRITE (6, 5) (X(K,J) , K=1 , 5 ) 

WRITE (6, 5) (X(K,J) ,K=6,9) 

WRITE (2,*) (X(K,J) , K=1 , 9 ) 

2000  CONTINUE 

C  ****  PLOT  FULL  DELTA  VS.  REF.  DELTA 

XLAB  =  '  COMPARE  FULL  DELTA  AND  REF.  DELTA' 

CALL  QPICTR (X, 18, NPTS, QY ( 3 , 11 ) ,QX( 1 ) ,QMOVE( 00 ) ,QXLAB( XLAB) , QLABE 


C  ***  PLOT  FULL  DELTA  vs.  FULL  STEP 

XLAB  =  '  COMPARE  FULL  DELTA  AND  FULL  STEP' 

CALL  QPICTR ( X , 1 8 , NPTS , QY ( 3 , 7 ) ,QX( 1 ) ,QMOVE( 00 ) , QXLAB ( XLAB ) , QLABEL 
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c  ****  PLOT  REF.  DELTA  vs.  ref  step 

XLAB  -  '  COMPARE  REF.  DELTA  AND  REF.  STEP' 

CALL  QPICTR(X, 18 ,NPTS,QY( 11 , 15) ,QX( 1 ) ,QMOVE( 00) , QXLAB ( XLAB ) , QLAB 


C  ****  PLOT  FULL  STEP  VS.  REF  STEP 

XLAB  =  '  COMPARE  FULL  STEP  AND  REF.  STEP' 

CALL  QPICTR(X, 18,NPTS ,QY( 7, 15 ) ,QX( 1 ) ,QMOVE( 00) , QXLAB (XLAB) , QLABE 


C  *****  PLOT  REAL  FREQUENCIES  FOR  FULL  DELTA  AND  STEP  FUNCTION  EQUILBRIA 
XLAB  •=  '  COMPARE  REAL  FREQ.  FOR  DELTA  &  STEP' 

CALL  QPICTR( X ,  18  , NPTS ,QY(5,9,13,17) ,QX(1) ,QMOVE(00) , QXLAB (XLAB) ) 


5  FORMAT ( 5D14 . 4 ) 

STOP 

END 

SUBROUTINE  FNl(Z,FZ) 

DOUBLE  COMPLEX  Z , FZ ( 2 ) , DLD , DTPlD , DTMlD , CHIBD , CHIAD , OMEG 
DOUBLE  COMPLEX  ZF,ZFP 

DOUBLE  PRECISION  DENS , WIG , LAM , BETA, BETAl , PERP , PERP1 , GAM 
DOUBLE  PRECISION  ALPl , ALP2 , DEL 

COMMON  DENS , WIG, LAM, BETA, BETAl , PERP, PERPl , GAM, ALPl , ALP2 , DEL 
OMEG  =  (Z  -BETA* LAM) /ALPl 
CALL  Z ETA ( OMEG , Z F , Z FP ) 

ZF  «  OMEG*OMEG*ZFP/( DEL* DEL* ( 1 . —BETA* BETA ) ) 

DLD  «  LAM* LAM* OMEG* OMEG  -  DENS*ZF 

DTPlD  -  Z*Z  - ( LAM+1 . ) * ( LAM+1 . ) -DENS* ( 1 . -PERP/2 . ) 

DTMlD  =  Z*Z  -(LAM-1. )*( LAM-1. )-DENS*(l. -PERP/2. ) 

CHIBD  =  DENS* ( BETA* DEL* OMEG* (1.-3. * PERP/2 . ) + PERP/2 . -1 . ) *ZF 
CHIAD  =  DENS*OMEG*OMEG*( 1 .-PERP/2. )  -  DENS* (( 1 . -PERP/2 .) * 

1  (1. -PERP/2. )-2.*BETA*DEL*OMEG*(l. -PERP/2. ) * ( 1 . -3*PERP/2 . ) )*ZF 
FZ ( 1 )  =DLD* DTPlD* DTMlD  +  ( WIG/2 .)*( DTPlD  +DTM1D ) * ( LAM*LAM*CHIAD 
1  -DENS* DENS* ( DEL* DEL* BETA* BETA* (1.-3. * PERP/2 . ) * ( 1 . -3 . *PERP/2 . ) 
1  *ZF*ZF  +(1.-  PERP/2 .) *ZF) ) 

FZ  (  2 )  =  DENS *ZF*DLD*DTMlD/ ( LAM* LAM )  -  (WIG/2 .) *CHIBD*CHIBD 

RETURN 

END 

SUBROUTINE  FN2(Z,FZ) 

DOUBLE  COMPLEX  Z , FZ ( 2 ) , DLS , DTPl S , DTMlS , CHIBS , CHI AS  ,  OMEG 
DOUBLE  COMPLEX  ZF,ZFP 

DOUBLE  PRECISION  DENS , WIG , LAM , BETA , BETAl , PERP , PERPl ,  GAM 
DOUBLE  PRECISION  ALPl , ALP2 , DEL 

COMMON  DENS , WIG, LAM , BETA, BETAl , PERP , PERPl , GAM , ALPl , ALP2 , DEL 
OMEG  =  (Z  -BETAl*LAM* ( 1 . -PERPl/4 . ) )/( DEL*GAM*LAM*ALP2 ) 

CALL  ZETA ( OMEG , ZF , ZFP ) 

ZF  =  OMEG*OMEG*ZFP/(DEL*DEL*ALP2 ) 

DLS  =  LAM * LAM *OMEG* OMEG  -  DENS*ZF/GAM 

DTPlS=  Z*Z-( LAM+1. )*( LAM+1. ) -DENS* GAM* ( 1 . -PERP 1/2 .+PERPl*PERPl/4 

DTMl S=  Z*Z-( LAM-1 . )*( LAM-1.  ) -DENS  * GAM* ( 1 . -PERP 1/2 . + PERPl * PERP 1/4 

CHIBS  =DENS* ( BETAl *DEL*GAM*OMEG* (1.-3. *PERPl/2 . ) +PERP1/2 . -1 . )*ZF 
CHI AS  =DENS*GAM*GAM*GAM* (1.-9. *PERPl/4 . ) *OMEG*OMEG  -DENS*GAM* ( 

1  ( 1 . -PERPl/2 . ) *( 1 .-PERPl/2 . ) -2 . * BETAl *DEL*GAM*OMEG* ( 1 . -PERP 1/2 . 

1  *( 1 .-3 . *PERPl/2 . ) ) *ZF 

FZ  (  1 )  =  DLS *DTPlS*DTMl S+ ( WIG/2 .  ) * ( DTPl S  +  DTMl S ) * ( LAM* LAM* CHI AS 
1  - DENS * DENS * GAM *GAM* ( BETAl * BETAl * DEL* DEL* ( 1 . -3 . *PERPl/2 . )* 

1  ( 1 . -3 . *PERPl/2 . ) *ZF*ZF  +  ( 1 . -9 . *PERPl/4 . ) *ZF ) ) 

FZ ( 2 )  =  DENS*ZF*DLS*DTMlS/( LAM* LAM* GAM )  - (WIG/2 .) *CHIBS*CHIBS 
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RETURN 

END 

SUBROUTINE  MULLER ( FN , FNREAL , ZEROS , N , NPREV, MAXI T , EPl , EP2 , M ) 

C  DETERMINES  UP  TO  N  ZEROS  OF  THE  FUNCTION  SPECIFIED  BY  FN  USING 
C  QUADRATIC  INTERPOLATION,  i.e.  MULLER'S  MEHTOD 
EXTERNAL  FNl , FN2 
LOGICAL  FNREAL 

INTEGER  MAXIT , N , NPREV , KOUNT , L , M 
DOUBLE  PRECISION  EPl , EP2 , EPS1 , EPS2 

DOUBLE  COMPLEX  ZEROS{N) , C , DEN, DIVDFl , DIVDF2 , DVDF1P , FZR( 2 ) , FZRDFL 
1 , FZRPRV,H, ZERO, SQR,HPREV, FN 

C  *********************  INPUT  ***************************** 

C  FN  NAME  OF  SUBROUTINE,  OF  THE  FORM  FN(X,FX)  WHICH  FOR  GIVEN  X 
C  RETURNS  F ( X ) ,  THIS  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT  IN  MAIN 
C  CALLING  PROGRAM. 

C  FNREAL  IS  A  LOGICAL  VARIABLE,  IF  .TRUE.  ALL  APPROX.  ARE  TAKEN 
C  TO  BE  REAL  ,  ALLOWING  THIS  ROUTINE  TO  BE  USED  EVEN  IF  F(X)  IS  ONLY 
C  DEFINED  FOR  REAL  X. 

C  ZEORS ( 1 ) . . . .  ZEROS (NPREV)  CONTAINS  PREVIOUSLY  FOUND  ZEROS  OF  THE 
C  FUNCTION,  PROVIDED  NPREV  .GT.  0 

C  ZEROS ( NPREV+1 )... .  ZEROS(N)  CONTAINS  FIRS  GUESS  FOR  THE  ZEROS 
C  TO  BE  FOUND 

C  MAXIT  IS  THE  MAXIMUM  NUMBER  OF  FUNCTION  EVALUATIONS  ALLOWED  @  ZERO. 
C  EPl  ITERATION  IS  STOPPED  IF  ABS(H)  . LT .  EPl*ABS(ZR),  WITH 
C  H  EQUAL  TO  THE  LATEST  CHANGE  IN  THE  ZERO  ESTIMATE 

C  EP2  ALTHOUT  THE  EPl  CRITERION  IS  NOT  MET,  ITERATION  IS  STOPPED  IF 
C  ABS ( F ( ZERO ) )  .LT.  EP2 

C  N  IS  THE  TOTAL  NUMBER  OF  ZEROS  TO  BE  FOUND 
C  NPREV  IS  THE  NUMBER  OF  ZEROS  FOUND  PREVIOUSLY 
C********************  OUTPUT  ************************** 

C  ZEROS ( NPREV  +1)  ...  ZEROS(N)  APPROXIMATIONS  TO  ZEROS 

C  INITIALIZTION 

EPSl  =  DMAXl (EPl,l.D-12) 

EPS2  =  DMAXl ( EP2, l.D-20) 

c 

DO  500  I =NPREV  +1,N 
KOUNT  =  0 

C  COMPUTE  FIRST  THREE  ESTIMATES  FOR  ZERO  AS  ... 

C  ZEROS ( I ) + . 5 ,  ZEROS ( I )  -  .5,  ZEROS(I) 

401  ZERO  =  ZEROS ( I ) 

H  =  .  5 

CALL  DFLATE ( FN , ZERO+ .5,1, KOUNT , FZR , DVDFlP , ZEROS , L , M ) 

I F ( L  .NE.  0)  GO  TO  401 

CALL  DFLATE ( FN , ZERO- .5,1, KOUNT , FZR , FZRPRV , ZEROS , L , M ) 

I F ( L  .NE.  0)  GO  TO  401 
HPREV  =  -1. 

DVDFlP  =  (FZRPRV  -  DVDF1 P ) /HPREV 

CALL  DFLATE ( FN , ZERO , I , KOUNT , FZR , FZRDFL , ZEROS , L , M ) 

I F ( L  .NE.  0)  GO  TO  401 

C  DO  WHILE  KOUNT  .LE.  MAXIT  OF  H  IS  RELATIVELY  BIG 
C  OR  FZR  =  F ( ZERO )  IS  NOT  SMALL 

C  OR  FZRDFL  =  FDFLATED( ZERO)  IS  NOT  SMALL  OR  NOT  MUCH 

C  BIGGER  THAN  ITS  PREVIOUS  VALUE  FZRPRV. 

440  DIVDFl  =  (FZRDFL  -  FZRPRV)/H 

DIVDF2  =  (DIVDFl  -  DVDFlP)/(H  +  HPREV) 

HPREV  =  H 

DVDFlP  =  DIVDFl 

C  =  DIVDFl  +  H*DI VDF2 

SQR  =  C*C  -  4 . * FZRDFL*DI VDF2 

IF  (FNREAL  .AND.  DREAL(SQR)  .LT.  0.)  SQR  =  (0.,0.) 

SQR  =  CDSQRT(SQR) 

IF  ( DREAL ( C ) *DREAL ( SQR ) +DIMAG ( C ) *DIMAG ( SQR )  .LT.  0.)  THEN 
DEN  =  C  -  SQR 
ELSE 

DEN  =  C  +  SQR 
END  IF 
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I F ( CDABS ( DEN )  .LE.  0.)  DEN  — (l.,0.) 

H  -  -2 . *FZRDFL/DEN 
FZRPRV  =  FZRDFL 
ZERO  **  ZERO  +  H 

I F ( KOUNT  .GT.  MAXIT)  GO  TO  499 
470  CALL  DFLATE(FN, ZERO, I, KOUNT, FZR, FZRDFL, ZEROS, L,M) 

I F( L  .NE.  0)  GO  TO  401 
C  CHECK  FOR  CONVERGENCE 

I F ( CDABS ( H )  .LT.  EPSl*CDABS ( ZERO ) )  GO  TO  499 
IF(  niIAXl(  CDABS  (  FZr.(M)  ),  CDABS{  FZRDFL)  )  .LT.  EPS2  )  GO  TO  49a 
C  CHECK  FOR  DIVERGENCE 

IF( CDABS ( FZRDFL)  .GE.  1 0 . *CDABS ( FZRPRV ) )  THEN 
H  =  H/2 . 

ZERO  =  ZERO  -  H 
GO  TO  470 
ELSE 

GO  TO  440 
END  IF 

499  ZEROS ( I )  =  ZERO 

500  CONTINUE 
RETURN 
END 

SUBROUTINE  DFLATE{ FN , ZERO , I , KOUNT, FZERO , FZRDFL , ZEROS , L , M ) 

C  TO  BE  CALLED  BY  MULLER 

INTEGER  I , KOUNT, J,L,M 

DOUBLE  COMPLEX  FZERO( 2 ), FZRDFL , ZERO , ZEROS ( 8 ), DEN 
L=  0 

KOUNT  =  KOUNT  +  1 
CALL  FN( ZERO, FZERO) 

FZRDFL  =  FZERO(M) 

I F ( I  .LT.  2)  RETURN 

DO  410  J=2 , I 

DEN  =  ZERO  -  ZEROS ( J-l  ) 

IF( CDABS (DEN)  . EQ .  0.)  THEN 

ZEROS ( I )  «  ZERO  *  1.001 

L=  1 

RETURN 

ELSE 

FZRDFL  =  FZRDFL/DEN 
END  IF 

410  CONTINUE 

RETURN 
END 

SUBROUTINE  SORTl ( N , ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 

DOUBLE  COMPLEX  ZEROS (N) 

DOUBLE  PRECISION  DUMl , DUM2 , DUM3 , DUM4 
INTEGER  I , J , K , N 

C  FIND  THE  LARGEST  GROWTH  RATE  IN  THE  DATA  SET  ZEROS(N) 

DUMl  =  DIMAG( ZEROS ( 1 )  ) 

J  =  1 

DO  100  I  =  2 , N 

IF  {DUMl  .GT.  DIMAG( ZEROS ( I ) ) )  THEN 
DUMl  =  DUMl 
J  =  J 
ELSE 

DUMl  =  D I MAG ( ZEROS (I) ) 

J  =  I 
END  IF 

100  CONTINUE 

DUM2  =  DREAL ( ZEROS ( J ) ) 

I F ( DUMl  .LE.  0.)  DUMl  =  0. 

C  DUMl  IS  THE  LARGEST  GROWTH  RATE  IN  THE  DATA  SET 

C  DUM2  IS  THE  REAL  FREQUENCY  CORRESPONDING  TO  THE  MAX  GROWTH 

C  NEXT  FIND  THE  NEXT  LARGEST  GROWTH  RATE  THAT  IS  NOT  EQUAL  IN 
C  MAGNITUDE  TO  THE  FIRST  GROWTH  RATE 
I F ( J  .EQ.  1)  J  =  3 
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M  =  J  -  1 

DUM3  =  DIMAG ( ZEROS { M ) ) 

DO  300  1  =  1, N 

IF(  DIMAG( ZEROS ( I ) )  . GE .  DUM3  .AND.  DIMAG ( ZEROS ( I ) )  .NE. 

1  DUMl )  THEN 

DUM3  =  DIMAG ( ZEROS ( I ) ) 

K  =  I 
END  IF 
CONTINUE 

DUM4  =  DREAL ( ZEROS ( K ) ) 

DUM3  IS  THE  SECOND  LARGEST  DISTINCT  GROWTH  RATE 

DUM4  IS  THE  REAL  FREQUENCY  CORRESPONDING  TO  THE  ABOVE 

RETURN 

END 

subroutine  zeta ( argl , arg2 , arg3 ) 

plasma  dispersion  funcion.  see  fried  and  conte  for  definition, 
received  from  d.g.  swanson,  8/12/64  cal  teck.  tested  by  martha 
pennell,rle,  mit.  shifted  into  fortran  iv  by  m.  liedberman.  1/67 
a.b.  lagdon,  6/69-  rewrote  extensively,  is  now  intelligible  and 
much  faster,  accuracy  same  for  small  arguments,  derivative  now 
works  for  large  arguments. 


argl  =  argument  of  z  function 

arg2  =  value  of  z  function 

arg3  =  value  of  derivative  of  z  function 


double  complex  argl , arg2 , arg3 

double  complex  z , al , a2 , a3 , bl , b2 , b3 , zl , zsq , aul , au5 
double  complex  bb 

double  precision  x , y , yabs , aa , daa , e r ror , ul , u2 , u3 , u4 , u5 

double  precision  sl,s2,fn 

equivalence  (bb,rbb) 

z  =  argl 

x  =  dreal(z) 

y  =  dimag(z) 

yabs  =  dabs(y) 

if(yabs  .It.  1.)  goto  10 


continued  fraction  method 

z  =  dcmplx ( x , yabs ) 

aa  =  0 . 

daa  =  1.5 

bb  =  1.5  -  z*z 

al  =  0. 

a2  =  -1. 

bl  =  1. 

b2  =  bb 

aul  =  a2/b2 

aa  =  aa  -  daa 

daa  =  daa  +  2. 

rbb  =  rbb  +  2. 

a3  =  a2*bb  +  al*aa 

b3  =  b2*bb  +  bl*aa 

zl  =  a3/b3 

au5  =  zl  -  aul 

i f ( dabs ( dreal ( au5 )) +dabs ( dimag ( au5 )). It .  l.d-7)  goto  5 

al  =  a2 

bl  =  b2 

a2  =  a3 

b2  =  b3 

aul  =  zl 

n  =  n  +  1 

goto  3 
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i f ( y )  52,51,51 
z  =  dcmplx(x,y) 
if(y*y-x*x  .ge.  -85.0)  then 

zl  =  dconjg(zl)  -( 0 ., 7 . 08981540362206 ) *z*cdexp( -z*z ) 
else 

zl  =  dconjg(zl) 
end  i  f 
arg3  =  zl 

arg2  =  -(1  +.5*zl)/z 
return 

if(abs(x)  .It.  4.)  goto  20 


asymptotic  series  method 
zsq  =  z*z 

zl  =  (0. ,1.77245385090551 ) * cdexp ( -z sq ) 
n  =  1 

aul  =  0.5/zsq 
au5  =  1.0/z 

au5  =  au5*aul *df lot j ( n ) 
zl  =  zl  -  au5 

error  =  dabs ( dreal ( au5 )) +dabs ( dimag ( au5 ) ) 
n  =  n  +  2 

if(error  .gt.  l.d-7)  goto  11 
arg3  =  -2.*z*zl 
arg2  =  zl  -  1.0/z 
return 


power  series  method,  needs  double  precision  u's  and  s's  on  some 
computers,  change  abs  to  dabs,  real  to  dreal  aimag  to  dimag  also 
ul  =  -2.0*(x*x  -  y*y) 
u2  =  -4.0*x*y 

error  =  1 . d-7/( dabs ( dreal ( z )) +dabs ( dimag ( z )) +1 . d-7 ) 
si  =  1.0 
s2  =  0.0 
n  =  3 

u5  =  dflotj(n) 
u3  =  ul/u5 
u4  =  u2/u5 
si  =  si  +  u3 
s2  =  s2  +  u4 

i f ( dabs ( u3 ) +dabs ( u4 )  .It.  error)  goto  25 

n  =  n  +  2 

fn  =  dflotj(n) 

u5  =  ( u3*ul-u4*u2 )/f n 

u4  =  ( u4*ul+u3*u2 )/f n 

u3  =  u5 

go  to  21 

zl  =  ( 0 ., 1 . 77245385090551 ) *cdexp ( -z*z ) -2 . *z*dcmplx ( si , s2 ) 

arg3  =  -2. *(1 ,+zl*z) 

arg2  =  zl 

return 

end 
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TEST  PROGRAM  FOR  MODEL  EQUATIONS  IN  HIGH  POWER 
FEL  PROBLEM. 

REVISION  12/3  TO  CORRECT  BOUNDARY  CONDITIONS 

REVISION  1/14  TO  INCLUDE  MULTIPLE  MODES 

REVISION  1/15  TO  INCLUDE  GROWTH  OF  TOTAL  WAVE  ENERGY 

REVISION  1/16  TO  INCLUDE  PHENOMENOLOGICAL  DAMPING 

REAL  A( 3000,50) ,PSIM(4000,50) ,PSI0(4000,50) , EPS , BETAl , TAU , DUM ( 50 

REAL  TIME ( 3000 ) ,UOLD( 4000 , 50 ) ,UNEW( 4000,50) , GAM( 3000 ) 

REAL  GROV7TH(  4090 , 50  )  ,  KWIGR  ,  KWIGL  ,  NU  ,  RISE  ,  DUMl  (  50  )  ,OMEGA(50)  ,KPOD 

REAL  ARPRIM( 50) ,AIMPRIM( 50 ) , ARN ( 50 ) ,ARNl( 50 ) ,AIMN( 50 ) ,AIMNl( 50 ) 
REAL  GROE( 3000 ) ,EWAV( 3000 ) 

INTEGER  J , NPART , NTIMES , N , F , NMODE , NSEP 
PARAMETER  (PI-3.141592653589) 

CHARACTER* 40  GLAB , XLAB , YLAB 
CHARACTER* 80  ABl , AB2 ( 50 ) , AB3 , AB4 
PS00 ( M , J )  =  -OMEGA( M) * FLOAT ( J ) *TAU 
FORMAT (  'WAVE  AMPLITUDE  MULTI  W/  R=',F5.3) 

XLAB  ='  TIME$ ' 

YLAB  =  '  WAVE  AMPLITUDE$ ' 

WRITE( 6,101) 

FORMAT ( '  INPUT  NO.  OF  PART . , NO .  OF  INTERAT . , GAM , BETAWIG , 

1  KWIGr ,BUDKER,NWIG, EPS, F, RISE, NPLUS, NMODE, NSEP, REF  ') 

READ ( 5 , * ) NPART, NTIMES , GAM0 , BETAW, KWIGR , NU , NWIG , EPS , F , RISE , NPLUS 
1  , NMODE, NSEP, REF 
WRIT£( 6,103) 

FORMAT ( '  INPUT  DATA : NPART , NTIMES , GAM , BETAW , KWIGR , NU , NWIG 
1, EPS, F, RISE, NPLUS, NMODE, NSEP, REF  IS: ' ) 

WRITE ( 6, *) NPART, NTIMES, GAM0, BETAW, KWIGR, NU, NWIG, EPS, F, RISE, NPLUS 

1  , NMODE, NSEP, REF 

KWIGL  =  2 . * FLOAT ( NWIG ) *PI 

BETA0  =  SQRT ( 1 .  -  1 ./( GAM0*GAM0 ) ) 

BETAZ0  -  SQRT(BETA0*BETA0  -  BETAW* BETAW/2 . ) 

NOPT  -  JIFIX(2.*FLOAT(NWIG)*BETAZO/(1.-BETAZO) )  +NPLUS 
OMEGA ( 1 )  =  ( FLOAT ( NOPT )*PI )/BETAZ0 
KPOD(l)  =  KWIGL  +  BETAZ0*OMEGA( 1 ) 

DO  50  J=2 , ( NMODE-1 )/2  +  1 

OMEGA ( J )  -  FLOAT( NOPT+( J-l ) *NSEP ) *PI/BETAZ0 
KPOD(J)  »  KWIGL  +BETAZ0*OMEGA( J) 

OMEGA( NMODE+2-J )  =FLOAT( NOPT- ( J-l ) *NSEP ) *PI/BETAZ0 
KPOD(NMODE+2-J)  =  KWIGL  +  BETAZ0*OMEGA( NMODE+2-J ) 

CONTINUE 
WRITE( AB4 , 1 ) REF 

WRITE (ABl , 104 ) NPART, NTIMES ,GAM0 , BETAW, NWIG 

FORMAT ( 'NPART-' ,I4,2X, 'NTIMES-' ,I4,2X, 'GAM-' ,F7.5,2X, 

1  'BETAWIG-' ,F8. 4, 2X, 'NWIG-' ,13) 

FORMAT ( ' N— ' , 1 4 , 3X , ' NU- ' ,E10.4,3X, ' KPOD- '  , El 0.4 
1 , 3X , 'OMEGA-' , E10 . 4 , 2X, ' NP- ' ,14) 

WRITE ( AB3, 106) KWIGL, EPS, F, RISE, KWIGR 

FORMAT ( 'KWIGL-' ,F8.4,2X, 'EPS-' , El  0 . 4 , 2X , ' F= ' , 1 3 , 2X , 

1 ' RISE— ' ,F8.4,2X, 'KWIGR-' , F8 . 4 ) 

BETAl  -  4 . *NU* BETAW* ( KWIGL ) *  *  2/ ( BETAZ0  *  BETA0  *GAM0 
1  * ( KWIGR*KWIGR ) ) 

ALPHAl  -  -BETAW/( 2 . *BETAZ0*BETAZ0 ) 

ALPHAR2  =  (1.  -  REF )/BETAZ0 

ALPHAl 2  --4 . *NU*KWIGL* * 2 * ( 1 . -BETAW* * 2/2 . )/( BETAZ0**2*KWIGR**2 
1  *GAM0) 

INITIALIZE  PHASE  AND  AMPLITUDE 
TIME ( 1 )  =  1. 

GROE(l)  =  0. 

TAU  =1 ./FLOAT (NPART- 1 ) 

TAUT  -  FLOAT ( F ) *TAU 
EWAV0  =  0. 

DO  60  Jl-1, NMODE 
ARPRIM(Jl)  =  0. 
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AIMPRIM(Jl)  =  0. 

AIMN(Jl)  =  0. 

ARN(Jl)  =  EPS 
A(  1  ,  J1 )  =  EPS 

EWAVO  =  EWAVO  + ( GAMO *KWIGR ) * *2  * ( ( KPOD ( J 1 ) -KWIGL ) *A ( 1 , Jl ) / 

1  KWIGL ) *  *2/( 4 . *NU* ( GAMO-1 .  ) ) 

GROWTH (l,Jl)  -  0. 

WRITE( AB2 ( Jl ) , 105) JIFIX{OMEGA( Jl ) *BETAZ0/PI ) ,NU,KPOD( Jl ) 

1  , OMEGA ( Jl ) , NPLUS 
DO  100  J=1 , NPART 

PSI0(J,J1)  =  PS00(Jl,J-l)  +  (KPOD(Jl)  -  OMEGA ( J 1 ) ) * ( NPART  -  J)*T 

UOLD(J,Jl)  =  KPOD(Jl)  -  OMEGA ( Jl ) 

CONTINUE 
CONTINUE 
WRITE( 6,102) 

FORMAT ( '  THE  WAVE  AMPLITUDES  ARE:  ') 

BEGIN  LOOP  FOR  TIME  INCREMENTS 

DO  1000  N  =2 , NTIMES 

TIME ( N )  =  1.  +  FLOAT( N-l ) *TAUT 

TRISE  =  (1.  -EXP (  ( 1 .-TIME(N) )/RISE  )) 

DO  80  Jl=l , NMODE 
DUM(Jl)  =  0. 

DUMl(Jl)  =  0. 

CONTINUE 

BEGIN  PART  II:  STEP  AMPLITUDES  AND  PHASES 
COMPLETE  SUM  FOR  AMPLITUDE  STEP 
DO  200  J=l, NPART 

BETAZ  =  BETAZ 0* ( UOLD( J ,  1 )  +  OMEGA ( 1 ) )/KPOD( 1 ) 

GAM ( J ) =  SQRT(  (l.-BETAZ*BETAZ)/(l./(GAM0*GAM0)+BETAW*BETAW/2. )  ) 

DUM2  =  0. 

DO  250  Jl=l , NMODE 

DUM2  =  DUM2  + ( KPOD ( 1 ) *KPCD (  Jl  ) -OMEGA ( Jl ) * BETAZ  0  * BETAZ  0  * 

1  ( UOLD ( J , 1 ) +OMEGA( 1 ) ) ) * ( ARN ( Jl ) *COS ( PS 1 0 ( J , Jl )  )  + 

1  AIMN (Jl )*SIN(PSI0(J ,J1 ) ) ) +BETAZ0* BETAZ 0* ( UOLD( J  ,  1  ) 

1  +OMEGA( 1 )  )*<ARPRIM( Jl ) *SIN( PS I  0 ( J, Jl )  ) -AIMPRIM ( Jl ) * 

1  COS ( PSI0 ( J , Jl ) ) ) 

DUM(Jl)  =  DUM(Jl)  +  COS { PS I 0 ( J , Jl ) ) *TAU*GAM ( J ) 

DUMl(Jl)  =  DUMl(Jl)  +  SIN(PSI0( J, Jl) )*TAU*GAM( J) 

CONTINUE 

PSIM( J  ,  1  )  =  PSI0(J,1)  +  TAUT*UOLD( J , 1 ) 

UNEW ( J , 1 )  =  UOLD(J.l)  +  TAUT*ALPHAl*GAM( J ) *GAM( J ) *DUM2 
DO  275  Jl=2, NMODE 

PSIM( J , Jl )  »  KPOD (Jl)*PSIM(J,l) /KPOD ( 1 ) 

1  + ( KPOD( Jl ) *OMEGA( 1 )/KPOD( 1 )  -OMEGA ( Jl )) *TIME ( N ) 

UNEW(  J  ,  Jl  )  «  KPOD  {  Jl )  *  (  UNEW (  J  ,  1 )  +OMEGA  (  1  )  )  /KPOD  ( 1  )  -OMEGA (  J 1  ) 

CONTINUE 

CONTINUE 

EWAV(N)  -  0. 

DO  280  Jl=l, NMODE 

ARN1 ( Jl ) =ARN( Jl ) +TAUT* ( ( ALPHAI 2 *AIMN ( Jl ) +BETA1 *DUM ( Jl ) )*TRISE/ 

1  OMEGA ( Jl ) -ALPHAR2  *ARN ( J 1 )  ) 

AIMNl ( Jl ) =AIMN ( Jl ) +TAUT* (  ( -ALPHA 1 2  *ARN ( Jl ) +BETA1 *DUM1 ( Jl )  ) *TRISE 
1  /OMEGA ( Jl ) -ALPHAR2  *AIMN ( Jl )  ) 

ARPRIM( Jl )  * ( ALPHAI 2  *AIMN ( J 1 ) +BETA1 *DUM ( J 1 )  ) *TRI SE/OMEGA ( J 1 ) 

1  -ALPHAR2  *ARN ( Jl ) 

AIMPRIM( Jl )  - (-ALPHAI 2 *ARN ( Jl )+BETAl*DUMl ( Jl ) ) *TRI SE/OMEGA ( Jl ) 

1  -ALPHAR2*AIMN( Jl ) 

A( N , Jl )  =  SQRT ( ARNl ( Jl ) *ARNl ( Jl )  +  AIMNl ( Jl ) *AIMNl ( Jl ) ) 

GROWTH ( N , Jl )  =  2. ' ( A( N , Jl ) -A( N-l , Jl ) )/( TAUT* ( A ( N, Jl ) +A( N-l , Jl ) ) ) 
EWAV(N)  =  EWAV  (  N )  (-  (  GAMO  *KWIGR  )  *  *  2  *  (  (  KPOD  (  Jl  )  -  KWIGL  )  *  A(  N  ,  J 1  )  / 

1  KWIGL) **2/( 4 . *NU* ( GAM0-1  .)  ) 

CONTINUE 

END  OF  PART  II  A) 

BEGIN  BOOKEEPING 
DO  400  J=1 , NPART- F 
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DO  375  Jl=l , NMODE 
PSI0(J,J1)  =  PSIM( J+F , Jl ) 

UOLD ( J , J 1 )  =  UNEW ( J+F , J 1 ) 

CONTINUE 
CONTINUE 
DO  500  J  =  F ,  1 , — 1 
DUM2  =  0. 

DO  475  Jl=l, NMODE 

DUM2  =  DUM2  + ( KPOD( J 1 ) *  KPOD ( 1 ) -OMEGA ( J 1 ) *BETAZ0  *BETAZ0 
1  *KPOD( 1 ) ) * ( ARN (Jl)*COS(PSOO(Jl, NPART- J+ ( N-l ) *  F )  )  + 

1  AIMN (Jl)*SIN(PS00(Jl, NPART- J+ ( N-l ) *  F )  ) )  + 

1  BETAZO*BETAZO*KPOD( 1 ) * ( ARPRIM( Jl ) * SIN ( PS00 ( Jl , NPART- J+ ( N-l ) *F) 

1  -AIMPRIM ( J 1 ) *COS (  PSO 0 ( J 1 , NPART- J+ ( N-l ) *  F )  ) ) 

CONTINUE 

UOLD( NPART+l-J , 1 )  =KPOD( 1 )-OMEGA( 1 )+( J-l ) *TAU*ALPHAl *DUM2 
PS 1 0 ( NPART+1- J ,1)=PS00(1, NPART- J+ (N-1)*F)+(J-1 ) *TAU*( KPOD( 1 )- 
1  OMEGA ( 1 )  +UOLD( NPART+1- J,l)  )/2. 

DO  495  Jl=2, NMODE 

UOLD( NPART+l-J , Jl )  =  KPOD ( J 1 )*( UOLD ( NPART+1- J , 1 ) +OMEGA ( 1 )) / 

1  KPOD ( 1 )  -  OMEGA ( Jl ) 

PS I 0( NPART+l-J, Jl )  =KPOD( Jl )*PSI0( NPART+1- J,1 )/KPOD( 1 ) 

1  + ( KPOD ( Jl ) * OMEGA ( 1 ) /KPOD ( 1 ) -OMEGA( J 1 ) )*TIME(N) 

CONTINUE 

CONTINUE 

END  OF  PART  II  B)  BOOKEEPING 
DO  600  Jl=l, NMODE 

GROE(N)  =  ( EWAV ( N ) -EWAV0 ) /( TAUT*  EWAV0 ) 

ARN(Jl)  =  ARNl(Jl) 

AIMN(  Jl  )  =  AIMNM  Jl  ) 

CONTINUE 
EWAV0  =  EWAV(N) 

REPEAT  PHASE  AND  AMPLITUDE  INCREMENT 

FOR  NEXT  TIME  STEP 

CONTINUE 

DO  2000  Jl=l, NMODE 
YLAB  =  '  WAVE  AMPLITUDE? ' 

GLAB  =  '  WAVE  AMPLITUDE  VS.  TIME  MULTI?' 

CALL  ANOTAT { %REF ( XLAB ) , %REF ( YLAB ) , 1,0,0, DSHL) 

CALL  PWRIT ( 500 , 80 , %REF( ABl ) ,70,1,0,0) 

CALL  PWRIT ( 500,48, %REF( AB2 (Jl)), 72, 1,0,0) 

CALL  PWRIT(  500 , 16 ,  %REF( AB3 ) , 7 5 , 1 , 0 , 0 ) 

CALL  EZXY (TIME,A(l,Jl) , NT I ME S , %REF ( AB4 ) ) 

YLAB  =  '  GROWTH  RATE  ?' 

GLAB  =  '  GROWTH  RATE  NORMALIZED  TO  TRANSIT  TIME?' 

CALL  ANOTAT ( %REF ( XLAB ) , %REF( YLAB) ,1,0, 0 , DSHL) 

CALL  EZXY ( TIME , GROWTH ( 1 , Jl ) , NTIMES , %REF ( GLAB )  ) 

CONTINUE 

GLAB  =  'FIELD  ENERGY  DENSITY  VS.  TIME?' 

YLAB  =  '  ENERGY  DENSITY?' 

CALL  ANOTAT ( %REF( XLAB ) , %REF( YLAB) ,1,0, 0 , DSHL) 

CALL  EZXY(TIME, EWAV, NTIMES, %REF(GLAB) ) 

GLAB  ='  RATE  OF  CHANGE  OF  ENERGY?' 

YLAB  =  'GROWTH  RATE?' 

CALL  ANOTAT ( %REF( XLAB ) , %REF( YLAB) ,1,0,0, DSHL ) 

CALL  EZXY(TIME, GROE, NTIMES, %REF(GLAB) ) 

STOP 

END 


QQQQQQQQQQ  6666666666666666666666666666666666666666666666666666  QQQQQQQQQQ 

IQQQQQQQQQ  Digital  Equipment  Corporation  -  VAX/VMS  Version  V4.6  QQQQQQQQQQ 

QQQQQQQQQ  6666666666666666666666666666666666666666666666666666  QQQQQQQQQQ 


M 

M 

AAA 

RRRR 

AAA 

BBBB 

L 

EEEEE 

MM 

MM 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M  M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

RRRR 

A 

A 

BBBB 

L 

EEEE 

M 

M 

AAAAA 

R  R 

AAAAA 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

BBBB 

LLLLL 

EEEEE 

PPPP 

AAA 

RRRR 

TTTTT 

EEEEE 

M 

M 

PPPP 

P  P 

A 

A 

R 

R 

T 

E 

MM 

MM 

P  P 

P  P 

A 

A 

R 

R 

T 

E 

M  M 

M 

P  P 

PPPP 

A 

A 

RRRR 

T 

EEEE 

M 

M 

PPPP 

P 

AAAAA 

R  R 

T 

E 

M 

M 

P 

P 

A 

A 

R  P 

T 

E 

M 

M 

P 

P 

A 

A 

R 

R 

T 

EEEEE 

M 

M 

P 

FFFFF 

000 

RRRR 

f  t 

4 

4 

F 

0 

O 

R 

R 

f  f 

4 

4 

F 

0 

O 

R 

R 

4 

4 

FFFF 

0 

0 

RRRR 

r  r 

44444 

F 

0 

0 

R 

R 

t  r 

4 

#  . 

F 

0 

O 

R 

R 

f 

4 

•  . 

F 

000 

R 

R 

t 

4 

I 


I 

I 

File  _VC$DRB1 : [MARABLE. FELJPARTEMP . FOR; 4  (3928,1,0),  last  revised  on 

I7-OCT-1984  15:18,  is  a  27  block  sequential  file  owned  by  UIC  [MARABLE).  The 
ecords  are  variable  length  with  a  fixed  control  size  of  2  bytes  and  implied 
(CR)  carriage  control.  The  longest  record  is  73  bytes. 


lob  PARTEMP  (2004)  queued  to  LN03QUE  on  21-MAR-1988  14:20  by  user  MARABLE,  UIC 
PmARABLE],  under  account  4790  at  priority  100,  started  on  printer  LTA8 :  on 
21-MAR-1988  14:20  from  queue  VCLN03B. 


|  QQQQQQQQQ 
QQQQQQQQQQ 
QQQQQQQQQQ 


6666666666666666666666666666666666666666666666666666 
Digital  Equipment  Corporation  -  VAX/VMS  Version  V4.6 
6666666666666666666666666666666666666666666666666666 


QQQQQQQQQQ 

QQQQQQQQQQ 

QQQQQQQQQQ 


I 

I 


/ 


100  C  MAIN  PROGRAM  FOR  COMPARISON  OF  DELTA  FUNCTION  AND  STEP  FUNCTION 
200  C  DISPERSIONS  USING  MULLER'S  METHOD  TO  DETERMINE  THE  COMPLEX  PART 
300  C  OF  THE  El GEN FREQUENCIES 

400  INTEGER  N , NPREV , MAXIT , I , M , L , NPTS 

500  DOUBLE  PRECISION  EFl , EP2 , DENS , WIG , LAM , BETA , PERP , BETAl , PERPl , GAM 

600  DOUBLE  PRECISION  DUMl , DUM2 , DUM3 , DUM4 , X ( 9 , 50 ) 

700  DOUBLE  PRECISION  TS , TD , LINC , LINC2 , LINC3 , Yl ( 8 ) , OMEG 

800  DOUBLE  PRECISION  ALP1 , ALP2 , DEL , GAMBAR , OMEG0 

900  CHARACTER* 4 0  XLAB 

1000  DOUBLE  COMPLEX  ZEROS ( 6 ) , ROOT, PREVRT( 20 ) 

1100  EXTERNAL  FN1 , FN2 

1200  LOGICAL  FNREAL 

1300  COMMON  DENS, WIG, LAM, BETA, BETAl , PERP, PERPl , GAM, ALPl ,ALP2 

1400  OPEN  (UNIT  =1,  NAME= ' OUTFI LE '  , STATUS  ='NEW') 

1450  OPEN  (UNIT  =2,  NAME= ’ DATINI ’  , STATUS  ='NEW') 

1500  FNREAL=  .FALSE. 

1600  MAXIT  =  1000 

1700  EPl  =  l.D-13 

1800  EP2  =  l.D-13 

1900  WRITE (6,1) 

2000  WRITE (1,1) 

2100  1  FORMAT ( '  INPUT  FOLLOWING  DATA:  DENS , WIG , LAM , LINC , LINC2 , LINC3 , BET 
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1  , PERP, DEL, NPTS  ') 

READ ( 5 , * ) DENS , WIG , LAM , LINC , LINC2 , LINC3 , BETA , PERP , DEL , NPTS 
WRITE ( 1 , * )  DENS , WIG, LAM, LINC, LINC2 , LINC3 , BETA, PERP , DEL , NPTS 
WRITE ( 6,*)  DENS, WIG, LAM, LINC, LINC2,LINC3, BETA, PERP, DEL, NPTS 
GAMBAR  -  DSQRT( 1 .  +  BETA* BETA  +  PERP  +  WIG) 

BETA  =  BETA/GAMBAR 

PERP  =  PERP/( GAMBAR*GAMBAR ) 

WIG  =  WIG/ (GAMBAR* GAMBAR) 

DENS  =  DENS/GAMBAR 

DEL  =  DEL/GAMBAR 

BETAl  =  BETA/DSQRT( 1 . -PERP ) 

PERPl  =  PERP/( 1 .-PERP) 

GAM  =  1 ./DSQRT{ 1 .-PERP) 

TD  =  PERP*GAMBAR*( 1 .  -  .5*WIG) 

TS  =  . 5* PERPl * GAMBAR* ( 1 .  -  .5*WIG  -  PERPl/3 . ) /GAM 
WRITE (6,3)  TD , TS 
WRITE( 1,3)  TD , TS 

FORMAT ( '  DELTA  TEMP  =  ',D12.3,'  STEP  TEMP  =  ',Dl2.3) 

WRITE (6,2) 

WRITE ( 1,2) 

FORMAT ( '  LAM  2nd  ROOT  DLD  DTPlD  DTMlD  DLS 

1  DTP1S  DTMlS ' ) 

DO  1000  L  =  1 , NPTS 
I F (  L  .LE.  20)  LAM  =  LAM  +  LINC 

I F (  L  .GT.  20  .AND.  L  .LE.  30)  LAM  =  LAM  +  LINC2 
I F (  L  .GT.  30)  LAM  =  LAM  +  LINC3 
ALPl  »  DEL*LAM* ( 1 . -  BETA* * 2 ) 

ALP2  =  DEL*LAM*GAM*( ( 1 . -BETAl *  * 2 ) -PERPl * ( 1 . - 3 . *BETAl *  *  2 ) /4 .  ) 

X ( 1 , L )  =  LAM 
N  =  6 
NPREV  =  0 

ZEROS ( 1 )  =  DCMPLX ( 0 . , 0 . ) 

IF  (L  .GT.  1)  THEN 
ZEROS ( 1 )  =  PREVRT(l) 

ZEROS ( 2 )  =  PREVRT ( 2 ) 

ZEROS ( 3 )  =  PREVRT ( 3 ) 

ZEROS ( 4 )  =  PREVRT ( 4 ) 

ZEROS ( 5 )  =  PREVRT ( 5 ) 

ZEROS ( 6 )  =  PREVRT ( 6 ) 

END  IF 

CALL  MULLER( FNl , FNREAL , Z EROS , N , NPREV, MAXI T , EPl , EP2 , 1 ) 

FIND  EIGENVALUES  WITH  THE  LARGEST  GROWTH  RATES 
CALL  SORTl ( N, ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 
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Yl ( 1 )  -  DUM3 
Yl ( 2 )  =  DUM4 
X( 2 , L )  =  DUMl 
X( 3 , L )  -  DUM2 
PREVRT(l)  -  ZEROS ( 1 ) 

PREVRT(2)  «  ZEROS ( 2 ) 

PREVRT ( 3 )  -  ZEROS ( 3 ) 

PREVRT ( 4 )  =  ZEROS ( 4 ) 

PREVRT ( 5 )  =  ZEROS ( 5 ) 

PREVRT ( 6 )  =  ZEROS (  6  ) 

C  ******  X2  IS  THE  LARGEST  GROWTH  RATE  IN  THEN  DATA  SET  FOR  FULL  DELTA 
C  ******  X3  IS  THE  REAL  FREQUENCY  CORESPONDING  TO  THE  ABOVE 
C  *******  yl  js  THE  SECOND  LARGEST  DISTINCT  IMAGINARY  FREQUENCY 
C  *******  y2  IS  THE  REAL  FREQUENC  CORRESPONDING  TO  THE  ABOVE 
OMEG  =  DUM2  -  BETA*LAM 

Yl  (  3  )  =  LAM* *2* ( OMEG*  *  2  -  DENS* ( 1 . -BETA* *2 ) ) 

Yl ( 4 )  =  DUM2* * 2  -(LAM  +1.)**2  -  DENS* ( 1 . -PERP/2 . ) 

Yl ( 5 )  -  DUM2*  *2  - ( LAM- 1 . ) *  *  2  -  DENS*(1.  -  PERP/2.) 

NPREV  =  0 

ZEROS ( 1 )  =  DCMPLX (  0  .  ,  0  .  ) 

IF  (L  .GT.  1)  THEN 
ZEROS ( 1 )  =  PREVRT ( 7 ) 

ZEROS ( 2 )  =  PREVRT ( 8 ) 

ZEROS ( 3 )  =  PREVRT ( 9 ) 

ZEROS ( 4 )  -  PREVRT ( 10 ) 

ZEROS ( 5 )  =  PREVRT(ll) 

ZEROS ( 6 )  -  PREVRT{ 12 ) 

END  IF 

CALL  MULLER ( FN2 , FNREAL , ZEROS , N , NPREV, MAXIT , EPl , EP2 , 1 ) 

CALL  SORT1 ( N , ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 

X( 4 , L )  =  DUMl 
X( 5 , L )  -  DUM2 
PREVRT ( 7 )  =  ZEROS ( 1 ) 

PREVRT ( 8 )  =  ZEROS ( 2 ) 

PREVRT ( 9 )  =  ZEROS ( 3 ) 

PREVRT (10)  =  ZEROS ( 4 ) 

PREVRT (11)  *  ZEROS ( 5 ) 

PREVRT (12)  =  ZEROS ( 6 ) 

C****  X4  IS  THE  LARGEST  GROWTH  RATE  IN  THE  DATA  SET  FOR  FULL  STEP 
C  ****  X5  IS  THE  REAL  FREQUENCY  CORESPONDING  TO  THE  ABOVE 
OMEG  =  DUM2  -  BETAl * LAM* ( 1 .  -  PERPl/4.) 

Yl ( 6 )  =  DENS  * GAM* ( ( 1 . -BETAl *  *  2 ) -PERPl * (1.-3. * BETAl *  *  2 ) /4 .  ) 

Yl ( 6 )  =  LAM**2* ( OMEG**2  -  Yl ( 6 ) ) 

Yl ( 7 )  =  DUM2**2- ( LAM+1 . ) **2-DENS* ( 1 . -PERPl/2 . + ( PERPl/2 . )**2)*GAM 
Yl  (  8 )  =  DUM2**2-( LAM-1 .  )* *2-DENS* ( 1 . -PERPl/2 .  +  ( PERPl/2 .  )**2)*GAM 
N  «=  4 
NPREV  =  0 

ZEROS ( 1 )  -  DCMPLX( 0 . , 0 . ) 

IF  (L  .GT.  1)  THEN 
ZEROS ( 1 )  =  PREVRT (13) 

ZEROS ( 2 )  =  PREVRT (14) 

ZEROS ( 3 )  =  PREVRT ( 1 5 ) 

ZEROS ( 4 )  =  PREVRT (16) 

END  IF 

CALL  MULLER ( FNl , FNREAL , ZEROS , N , NPREV, MAXIT , EPl , EP2 , 2 ) 

CALL  SORTl(N, ZEROS, DUMl, DUM2 , DUM3 , DUM4 ) 

X( 6 , L)  =  DUMl 
X( 7 , L  )  =  BETA* LAM 
PREVRT ( 13)  =  ZEROS ( 1 ) 

PREVRT (14)  =  ZEROS ( 2 ) 

PREVRT (15)  =  ZEROS ( 3 ) 

PREVRT ( 16 )  =  ZEROS { 4 ) 

C  *****  X6  IS  THE  LARGEST  GROWTH  RATE  FOR  THE  REFERENCE  DELTA  FUNCTION 
C  ****  X7  IS  AN  APPROXIMATION  OF  THE  REAL  FREQUENCY  FOR  SMALL  DENSITIES 
NPREV  =  0 

ZEROS ( 1 )  =  DCMPLX ( 0  .  , 0  .  ) 
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IF  (L  .GT.  1)  THEN 
ZEROS  ( 1 )  **  PREVRT(  17  ) 

ZEROS( 2)  -  PREVRT( 18 ) 

ZEROS! 3)  -  PREVRT (19) 

ZEROS! 4)  -  PREVRT! 20) 

END  IF 

CALL  MULLER! FN2 , FNREAL , ZEROS , N , NPREV , MAXIT , EPl , EP2 , 2 ) 

CALL  SORTl ( N , ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 

X( 8 , L )  -  DUMl 

X( 9 , L )  -  BETAl*LAM 

PREVRT! 17)  =  ZEROS! 1) 

PREVRT! 18)  -  ZEROS! 2) 

PREVRT! 19)  =  ZEROS! 3) 

PREVRT! 20)  -  ZEROS! 4) 

C  ******  X8  IS  THE  LARGEST  GROWTH  RATE  FOR  THE  REFERENCE  STEP  FUNCTION 
C  *****  X9  IS  AN  APPROXIMATION  OF  THE  REAL  FREQUENCY  FOR  SMALL  DENSITIES 
WRITE (1,5)  LAM , (Yl(K) , K=1 , 8  ) 

WRITE! 6, 5)  LAM, ( Yl ( K ) , K=1 , 8 ) 

1000  CONTINUE 

WRITE! 1,6) 

WRITE! 6, 6) 

6  FORMAT! '  LAM  FULL  DELTA  FULL  STEP  REF  DELTA 

1  REF  STEP  ' ) 

DO  2000  J  ■*  1 , NPTS 
WRITE! 1,5) (X(K, J) ,K-1,9) 

WRITE (6,5)(X(K,J),K“1,9) 

WRITE (2,*)(X(K,J) , K=1 , 9 ) 

2000  CONTINUE 

C  ****  PLOT  FULL  DELTA  vs.  REF.  DELTA 

XLAB  =  '  COMPARE  FULL  DELTA  AND  REF.  DELTA' 

CALL  QPICTR ( X , 18 , NPTS , QY ( 3 , 11 ) ,QX( 1 ) ,QMOVE( 00 ) , QXLAB ( XLAB ) , QLABE 


C  ***  PLOT  FULL  DELTA  vs.  FULL  STEP 

XLAB  -  '  COMPARE  FULL  DELTA  AND  FULL  STEP' 

CALL  QPICTR ( X , 18 , NPTS , QY ( 3,7 ) ,QX( 1 ) ,QMOVE( 00 ) , QXLAB (XLAB) , QLABEL 


C  ****  PLOT  REF.  DELTA  vs.  REF  STEP 

XLAB  =  '  COMPARE  REF.  DELTA  AND  REF.  STEP' 

CALL  QPICTR (X, 18, NPTS, QY (11, 15) ,QX(1 ) ,QMOVE(00) , QXLAB (XLAB) , QLAB 


C  ****  PLOT  FULL  STEP  VS.  REF  STEP 

XLAB  =  '  COMPARE  FULL  STEP  AND  REF.  STEP' 

CALL  QPICTR! X, 18 ,NPTS , QY (7,15) ,QX(1) ,QMOVE(00) , QXLAB ( XLAB ) , QLABE 


C  *****  PLOT  REAL  FREQUENCIES  FOR  FULL  DELTA  AND  STEP  FUNCTION  EQUILIBRI 
XLAB  *  '  COMPARE  REAL  FREQ.  FOR  DELTA  &  STEP' 

CALL  QPICTR(X,18,NPTS,QY( 5,9 ,13, 17) ,QX( 1 ) ,QMOVE( 00) , QXLAB (XLAB) ) 


5  FORMAT ( 9D9 . 2 ) 

STOP 

END 

SUBROUTINE  FNl(Z,FZ) 

DOUBLE  COMPLEX  Z  ,  FZ ( 2  ) , DLD , DTP1D , DTMlD , CHI BD , CHI AD , OMEG 
DOUBLE  PRECISION  DENS , WIG , LAM , BETA , BETAl , PERP , PERPl , GAM 
DOUBLE  PRECISION  ALPl,ALP2 


18400  COMMON  DENS , WIG , LAM , BETA, BETAl , PERP , PERPl , GAM, ALPl , ALP2 

18500  OMEG  =  Z  -  BETA* LAM  +  ALPl *DCMPLX ( 0 . , 1 . ) 

19100  DLD  -  LAM*LAM* ( OMEG*OMEG-DENS* ( 1 . -BETA*BETA ) ) 

19200  DTPlD  »  Z*Z  -(LAM  +1.)*(LAM  +1.)  -DENS * ( 1 . -PERP/2 . ) 

19300  DTMlD  »  Z*Z  -(LAM  -l.)*(LAM  -1.)  -DENS* ( 1 . -PERP/2 . ) 

19400  CHI BD  =  DENS*LAM* ( BETA*OMEG* (1.-3. *PERP/2 . )-LAM* ( 1 . -BETA*BETA) 

19500  1*( 1 .-PERP/2. ) ) 

19600  CHI AD  =  DENS* ( 2 . * BETA* LAM* OMEG* (1. -PERP/2. )*(l.-3.*PERP/2. ) - 

19700  1  LAM*LAM* ( 1 . -BETA*BETA ) * ( 1 . -PERP/2 . ) * ( 1 . -PERP/2 . )  + 

19800  1  ( 1 . — 3 . *PERP ) * OMEG* OMEG ) 

19900  FZ ( 1 )  =DLD*DTPlD*DTMlD  +  ( WIG/2 .)*( DTPlD  +DTM1D ) * LAM* LAM* ( CHI AD 

20000  1 -DENS  * DENS  * ( (1.-3. *PERP ) * ( 1 . -BETA* BETA ) +BETA*BETA* (1.-3. * PERP/2 

.  ) 

20100  1* ( 1 . -3 . * PERP/2 . ) ) ) 

20200  FZ ( 2 ) =DENS* ( 1 . —BETA* BETA ) * DLD* DTMlD- ( WIG/2 . ) *CHIBD*CHIBD 

20300  RETURN 

20400  END 

20500  SUBROUTINE  FN2(Z,FZ) 

20600  DOUBLE  COMPLEX  Z , FZ ( 2 ), DLS , DTP1S , DTMlS , CHIBS , CHIAS , OMEG 

20700  DOUBLE  PRECISION  DENS , WIG , LAM, BETA, BETAl , PERP , PERPl , GAM 

20750  DOUBLE  PRECISION  ALPl,ALP2 

20800  COMMON  DENS , WIG , LAM , BETA, BETAl , PERP , PERPl , GAM , ALPl , ALP? 

20900  OMEG  =  Z  -  BETAl * LAM* (1.-PERP1/4. )  +  ALP2*DCMPLX( 0 . , 1 . ) 

21500  DLS=DENS*GAM*( ( 1 . -BETAl* BETAl ) -PERPl* (1.-3. * BETAl* BETAl ) /4 . ) 

21600  DLS  =  LAM*LAM*(OMEG*OMEG  -  DLS) 

21700  DTPlS  «Z*Z— (LAM+1. )*(LAM+1. ) -DENS* GAM* ( 1 . -PERP 1/2 .+PERPl*PERPl/4 

-  ) 

21800  DTMlS  =Z  *Z- ( LAM-1 . )*(LAM-1. ) -DENS* GAM* ( 1 . -PERP 1/2 .+PERPl*PERPl/4 

.  ) 

21900  CHIBS  =DENS*LAM*GAM*GAM* ( BETAl * OMEG* (1.-3. * PERP 1/2 . ) -LAM* 

22000  1  ( 1 . -PERPl/2 . ) * ( 1 . -BETAl * BETAl -PERPl * (1.-3. * BETAl * BETAl ) /4 . ) ) 

22200  CHIAS=DENS*GAM*GAM*GAM* ( 2 . * BETAl * LAM* OMEG* ( 1 . -PERPl/2 .  ) * ( 1  .  -3 . * 

22300  1  PERPl/2. ) -LAM* LAM* ( (1. -PERPl/2. )*( 1 . -PERPl/2 . )*( 1 . -BETAl *BETA1 

22400  1  PERPl* (1.-3. * BETAl * BETAl )/4 . ) )+( 1 .-9 . *PERPl/4 . ) *OMEG*OMEG ) 

22600  FZ ( 1 )  =  DLS*DTP1S*DTM1S+ ( WIG/2 . ) * ( DTP 1S+ DTMlS ) * LAM* LAM* ( CHIAS  - 

22700  1  DENS*DENS*GAM**4* ( ( 1 . -9 . *PERPl/4 . ) * ( 1 . -BETAl *BETAl-PERPl * ( 1 . - 3 


22800  1 *BETAl *BETAl )/4 . ) +BETA1 * BETAl * ( 1 . -3 . *PERPl/2 . ) * ( 1 . -3 . *PERPl/2 . ) 

)  ) 

23000  FZ ( 2 ) =DENS*GAM* ( 1 . -BETAl *BETAl -PERPl * ( 1.-3. *BETAl*BETAl )/4 . ) *DLS 

23100  l*DTMlS  -  (WIG/2. )*CHIBS*CHIBS 

23200  RETURN 

23300  END 

23400  SUBROUTINE  MULLER ( FN , FNREAL , ZEROS , N , NPREV , MAXI T , EPl , EP2 , M ) 

23500  C  DETERMINES  UP  TO  N  ZEROS  OF  THE  FUNCTION  SPECIFIED  BY  FN  USING 
23600  C  QUADRATIC  INTERPOLATION,  i.e.  MULLER'S  MEHTOD 
23700  EXTERNAL  FN1 , FN2 

23800  LOGICAL  FNREAL 

23900  INTEGER  MAXIT, N, NPREV, KOUNT,L,M 

24000  DOUBLE  PRECISION  EPl , EP2 , EPSl , EPS2 

24100  DOUBLE  COMPLEX  ZEROS ( N ) , C , DEN , DIVDFl , DI VDF2 , DVDFlP , FZR ( 2 ) , FZRDFL 

24200  1 , FZRPRV, H, ZERO, SQR,HPREV, FN 

24300  C  *********************  INPUT  ***************************** 

24400  C  FN  NAME  OF  SUBROUTINE,  OF  THE  FORM  FN(X,FX)  WHICH  FOR  GIVEN  X 
24500  C  RETURNS  F(X),  THIS  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT  IN  MAIN 

24600  C  CALLING  PROGRAM. 

24700  C  FNREAL  IS  A  LOGICAL  VARIABLE,  IF  .TRUE.  ALL  APPROX.  ARE  TAKEN 

24800  C  TO  BE  REAL  ,  ALLOWING  THIS  ROUTINE  TO  BE  USED  EVEN  IF  F(X)  IS  ONLY 

24900  C  DEFINED  FOR  REAL  X. 

25000  C  ZEORS(l)....  ZEROS(NPREV)  CONTAINS  PREVIOUSLY  FOUND  ZEROS  OF  THE 
25100  C  FUNCTION,  PROVIDED  NPREV  .GT.  0 

25200  C  ZEROS ( NPREV+1 )... .  ZEROS(N)  CONTAINS  FIRS  GUESS  FOR  THE  ZEROS 
25300  C  TO  BE  FOUND 

25400  C  MAXIT  IS  THE  MAXIMUM  NUMBER  OF  FUNCTION  EVALUATIONS  ALLOWED  @  ZERO. 
25500  C  EPl  ITERATION  IS  STOPPED  IF  ABS(H)  .LT.  EPl*ABS(ZR),  WITH 
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C  H  EQUAL  TO  THE  LATEST  CHANGE  IN  THE  ZERO  ESTIMATE 

C  EP2  ALTHOUT  THE  EPl  CRITERION  IS  NOT  MET,  ITERATION  IS  STOPPED  IF 
C  ABS ( F ( ZERO ) )  .LT.  EP2 

C  N  IS  THE  TOTAL  NUMBER  OF  ZEROS  TO  BE  FOUND 
C  NPREV  IS  THE  NUMBER  OF  ZEROS  FOUND  PREVIOUSLY 
C********************  OUTPUT  ************************** 

C  ZEROS (NPREV  +1)  ...  ZEROS(N)  APPROXIMATIONS  TO  ZEROS 

C  INITIAL I ZTION 

EPS1  -  DMAX1 (EPl,l.D-12) 

EPS2  =  DMAX1 ( EP2 , 1 . D-2  0 ) 
c 

DO  500  I=NPREV  +1,N 
KOUNT  =  0 

C  COMPUTE  FIRST  THREE  ESTIMATES  FOR  ZERO  AS  ... 

C  ZEROS ( I )  +  . 5 ,  ZEROS ( I  )  -  .5,  ZEROS(I) 

401  ZERO  =  ZEROS ( I ) 

H  =  .  5 

CALL  DFLATE ( FN , ZERO+ .5,1, KOUNT , FZR , DVDFlP , ZEROS , L , M ) 

I F ( L  .NE.  0)  GO  TO  401 

CALL  DFLATE ( FN , ZERO- .5,1, KOUNT , FZR , FZRPRV , ZEROS , L , M ) 

I F ( L  .NE.  0)  GO  TO  401 
HPREV  -  -1. 

DVDFlP  =  (FZRPRV  -  DVDFl P ) /HPREV 

CALL  DFLATE ( FN , ZERO , I , KOUNT , FZR , FZRDFL , ZEROS , L , M ) 

I F ( L  .NE.  0)  GO  TO  401 

C  DO  WHILE  KOUNT  .LE.  MAXIT  OF  H  IS  RELATIVELY  BIG 
C  OR  FZR  -  F ( ZERO )  IS  NOT  SMALL 

C  OR  FZRDFL  -  FDFLATED( ZERO )  IS  NOT  SMALL  OR  NOT  MUCH 

C  BIGGER  THAN  ITS  PREVIOUS  VALUE  FZRPRV. 

440  DIVDFl  =  (FZRDFL  -  FZRPRV)/H 

DIVDF2  -  (DIVDFl  -  DVDFlP )/(H  +  HPREV) 

HPREV  =  H 

DVDFlP  -  DIVDFl 

C  »  DIVDFl  +  H*DIVDF2 

SQR  -  C*C  -  4 . *FZRDFL*DIVDF2 

IF  (FNREAL  .AND.  DREAL(SQR)  . LT .  0.)  SQR  =  ( 0 . , 0 . ) 

SQR  -  CDSQRT(SQR) 

IF  ( DREAL ( C ) *DREAL ( SQR ) +DIMAG ( C ) *DIMAG ( SQR )  .LT.  0.)  THEN 
DEN  =  C  -  SQR 
ELSE 

DEN  =  C  +  SQR 
END  IF 

I F ( CDABS ( DEN )  . LE .  0.)  DEN  =(1.,0.) 

H  =  -2 . *FZRDFL/DEN 
FZRPRV  =  FZRDFL 
ZERO  =  ZERO  +  H 
I F ( KOUNT  .GT.  MAXIT)  GO  TO  499 
470  CALL  DFLATE(FN, ZERO, I, KOUNT, FZR, FZRDFL, ZEROS, L,M) 

I F ( L  .NE.  0)  GO  TO  401 
C  CHECK  FOR  CONVERGENCE 

I F ( CDABS ( H )  .LT.  E PS 1  * CDABS ( ZERO ) )  GO  TO  499 
IF( DMAXK CDABS ( FZR (M) ) , CDABS ( FZRDFL ) )  . LT .  EPS2 )  GO  TO  499 

C  CHECK  FOR  DIVERGENCE 

IF(CDABS( FZRDFL)  .GE.  1 0 . *CDABS ( FZRPRV ) )  THEN 
H  =  H/2 . 

ZERO  =  ZERO  -  H 
GO  TO  470 
ELSE 

GO  TO  440 
END  IF 

499  ZEROS ( I )  =  ZERO 

500  CONTINUE 
RETURN 
END 

SUBROUTINE  DFLATE ( FN , ZERO , I , KOUNT , FZERO , FZRDFL , ZEROS , L , M ) 

C  TO  BE  CALLED  BY  MULLER 
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INTEGER  I , KOUNT , J , L , M 

DOUBLE  COMPLEX  FZERO ( 2 ) , FZRDFL , ZERO , ZEROS ( 8 ) , DEN 
L=  0 

KOUNT  =  KOUNT  +  1 
CALL  FN( ZERO, FZERO) 

FZRDFL  -  FZERO(M) 

I F ( I  .LT.  2)  RETURN 

DO  410  J=2 , I 

DEN  =  ZERO  -  ZEROS ( J-l ) 

I F ( CDABS ( DEN )  .EQ.  0.)  THEN 

ZEROS ( I )  =  ZERO  *  1.001 

L=  1 

RETURN 

ELSE 

FZRDFL  =  FZRDFL/DEN 
END  IF 

410  CONTINUE 

RETURN 
END 

SUBROUTINE  SORTl ( N , ZEROS , DUMl , DUM2 , DUM3 , DUM4 ) 

DOUBLE  COMPLEX  ZEROS(N) 

DOUBLE  PRECISION  DUMl , DUM2 , DUM3 , DUM4 
INTEGER  I , J , K , N 

C  FIND  THE  LARGEST  GROWTH  RATE  IN  THE  DATA  SET  ZEROS(N) 

DUMl  =DIMAG( ZEROS { 1 ) ) 

J  =  1 

DO  100  I  =  2 ,  N 

IF  (DUMl  .GT.  DIMAG ( ZEROS ( I ) ) )  THEN 
DUMl  =  DUMl 
J  =  J 
ELSE 

DUMl  -  DIMAG( ZEROS ( I ) ) 

J  =  I 
END  IF 

100  CONTINUE 

DUM2  -  DREAL ( ZEROS ( J ) ) 

I F ( DUMl  .LE.  0.)  DUMl  =  0. 

C  DUMl  IS  THE  LARGEST  GROWTH  RATE  IN  THE  DATA  SET 

C  DUM2  IS  THE  REAL  FREQUENCY  CORRESPONDING  TO  THE  MAX  GROWTH 

C  NEXT  FIND  THE  NEXT  LARGEST  GROWTH  RATE  THAT  IS  NOT  EQUAL  IN 
C  MAGNITUDE  TO  THE  FIRST  GROWTH  RATE 
I F ( J  .EQ.  1)  J  =  3 
M  =  J  -  1 

DUM3  =  DIMAG( ZEROS ( M ) ) 

DO  300  1=1, N 

IF(  DIMAG ( ZEROS ( I ) ) . GE . DUM3 . AND . DIMAG ( ZEROS ( I )  ) .NE. 

1  DUMl)  THEN 

DUM3  =  DIMAG( ZEROS ( I ) ) 

K  =  I 
END  IF 

300  CONTINUE 

DUM4  =  DREAL ( ZEROS ( K) ) 

C  DUM3  IS  THE  SECOND  LARGEST  DISTINCT  GROWTH  RATE 

C  DUM4  IS  THE  REAL  FREQUENCY  CORRESPONDING  TO  THE  ABOVE 

RETURN 
END 
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RESTART  CODE 

CODE  TO  EVALUATE  TEMPORAL  EVOLUTION  OF  THE  SPRECTRA 
OF  UNSTABLE  MODES  IN  A  HELICAL  WIGGLER  FREE  ELECTRON  LASER 
DELETION  OF  FIRST  TRANSIT  TIME 

FIELD  AND  PARTICLE  EQUATIONS  ARE  EVOLVED  BY  ADAMS-BASHFORTH 
METHOD  WITH  INITALIZATION  BY  RUNGE-KUTTA  METHOD 
REFORMULATION  OF  THE  PARTICLE  PHASE  3/13 
CONVERSION  TO  CRAY  FORTRAN 

INCLUSION  OF  FREQUENCY  SHIFT  ERROR  CHECK  IN  ADAMS-BASHFORTH 
EQUATION  SOLVER 

REPLACE  EXPRESSION  FOR  THE  DERIVATIVES  WITH  THE  FUNCTIONAL 

EVALUATION  OF  THE  DIFFERENTIAL  EQUATION 

PLOT  DATA  AFTER  EVERY  10  CALCULATIONS 

OUTPUT  DATA  IF  TIME  LIMIT  IS  APPROACHED 

MODIFICATIONS  TO  PRODUCE  RESTART  DATA  8/1 

REAL  BETAZO , BETAO , KPOD( 20 ) ,OMEGA(20) , BETAZ , GAM 

REAL  CTHET( 20) , PS 1(20,3 000,5), U0 (3000, 5) ,TIM 

REAL  TEMP ( 3000 ) , TEMPI ( 3000 ) , ATEMP ( 20 ) , ATEMPl ( 20 ) , THETA ( 20,5) 
REAL  TP( 20 , 5 ) , TIME ( 5000 ) ,TPLOT(6) 

REAL  FREQ( 5000 , 20 ) ,EWAV(5000) ,A(20,5) ,APP(20,5) 

REAL  PHI1( 20, 5 ) , PHI 2 ( 20, 5) , KWIGL , PS  1 1 ( 600 , 6  ) ,U0I ( 600,6 ) 

, CORA( 20) , TPC ( 20) , PLAI (5000,20) 

REAL  KWIGR , NU , NUR , NUI , FI LL , PLAR ( 5000,20) 

REAL  Fll(20) , F12 ( 20 ) , F13 ( 20 ) , Fl4 ( 20 ) , F15 ( 20 ) , F21 ( 20) ,F22(20) 
REAL  F23(20) ,F24(20) ,F25(20) 

INTEGER  F , MM , JP , J , Jl , K , MAXIT, TLIM , NCOUNT, NLAST 
COMMON/BLKl/BETAZO ,GAM0 , EETAW, KPOD , OMEGA 
COMMON/BLK 3/TIM , NPART , N , RISE , TAU , BETAl , BETA2 , F , NUI , NUR 
COMMON/BLK4/  ALPHAl , ALPHA2 , NMODE 
COMMON/BLK 5/PS I ,U0,PHIl,PHl2 , A, THETA, APP , TP 
PARAMETER  (PI-3.1415926535) 

PS00(J1,J)  =  -OMEGA( Jl ) *FLOAT( J-l ) *TAU 
TRISE(N)  =  1.  -RISE* ( EXP( ( -TIM+1 . )/RISE )  -1.) 

OPEN( UNIT-1 , FILE- ' FILEl ' , FORM= ' UNFORMATTED ' , STATUS* ' OLD ' ) 

OPEN ( UNIT— 2 , FI LE= ' FI LE2 ' , FORM* ' UNFORMATTED ' , STATUS- ' OLD ' ) 

OPEN ( UNIT-3,  FILE-' FI LE3' , FORM- ' UNFORMATTED ' , STATUS- ' OLD ' ) 

OPEN (UNIT-4,  FILE-' FI LE4 ' , FORM- ' UNFORMATTED ' , STATUS- ' NEW ' ) 
OPEN(UNIT=5,FILE='FILE5' , FORM- 'UNFORMATTED' , STATUS- ' NEW' ) 

OPEN ( UN IT- 7 , FILE— ' FI LE7 ' , FORM- ' UNFORMATTED' , STATUS- ' NEW ' ) 

READ ( 1 )  TIME, FREQ, EWAV, PLAR, PLAI , KPOD, OMEGA, KWIGR, NU, 

GAM0 , BETAW , RISE, FILL, REF, EPS , PHASE , ERROR , BETAZ  0 , ERROR2 , 

PSI I ,U0I ,TPLOT, TLIM, BETAO 

READ ( 2 )  NWIG , NPART , F , NMODE , NPLUS , MAXIT , NSEP , NTIMES , NCOUNT , 

NLAST 

READ ( 3 )  PSI ,U0 , Fll , F12 , F13 , F21 , F22 , F2 3 , A , APP , THETA, 

TP , PHI 1 , PHI 2 
NCOUNT  -  NCOUNT  +  1 

FORMAT (  '  THE  PRED.-CORR.  METHD  FAILED  TO  CONVERGE  ON  STEP',2X, 

17,'  AFTER  ' , I  4 , 2X, '  INTERATIONS ' , 2X , ' ON  MODE  ',13) 

KWIGL  =  2 . * FLOAT ( NWIG) *PI 

BETAl  -  2 . *FILL*NU*KWIGL**2*BETA0*BETAW/( KWIGR**2*BETAZ0**3 ) 

NUR  =  ( 1 .-REF) /BETAZO 

NUI  =  -4 . *FILL*NU*KWIGL**2* ( 1 . -BETAW* *2/2 . )/( GAM0* BETAZO* 
KWIGR**2*BETAZO*OMEGA( 1) ) 

BETA2  =  8 . *NU*BETA0*KWIGL**2/( BETAZ0*KWIGR* *  2 ) 

BETA2  =  0. 

ALPHAl  =  KPOD( 1 )/BETAZ0 

ALPHA2  =  . 5*BETAW*GAM0 * KPOD ( 1 )/BETAZ0**2 
TAU  =  1 ./FLOAT (NPART  -1) 

TAUT  =  FLOAT ( F ) *TAU 

NOW  EVOLVE  PARTICLES  AND  FIELDS  WITH  ADAMS-BASHFORTH  PREDICTOR 
CORRECTOR  METHOD  USING  THE  RESULTS  OF  THE  FOUR  PREVIOUS  TIMES 
AS  INITIAL  CONDITIONS 
DO  7000  N-NLAST+1 , NTIMES 
TIM  =  1.  +  FLOAT ( N-l ) *TAUT 
DO  5100  Jl-1, NMODE 
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CALL  EVOLV( Jl , F24 , F14 , 2 ) 

DO  5200  JP=1 ,NPART-4*F 
CALL  F4 ( JP+F, FOUT, 2 ) 

CALL  F4 ( JP+2*F , FOUTl , 3 ) 

CALL  F4( JP+3*F, FOUT2, 4 ) 

CALL  F4(JP+4*F, FOUT 3 , 5 ) 

U0<JP,1)  =  U0 ( JP+F , 2 )  +  TAUT* ( 55. * FOUT  -59.*FOUTl  +37.*FOUT2 
-9 . * FOUT 3 )/24 . 

PSI(1,JP,1)  =  PSI ( 1 , JP+F, 2 )  +  TAUT* ( 55.*U0(JP+F,2)  - 
59.*U0( JP+2*F, 3)+37.*U0( JP+ 3  * F , 4 ) -9 . *U0 ( JP+4*F,5) )/24. 

TEMP ( JP )  =  19 . *  FOUT  -5.*FOUTl  +  FOUT2 

TEMPI ( JP )  =  19.*U0(JP+F,2)-5.*U0(JP+2*F,3)+U0(JP+3*F,4) 

DO  5300  JP=NPART-4*F+1 , NPART-F 
CALL  F4( JP+F, FOUT, 2) 

U0(JP,1)  =  U0 { JP+F , 2 )  +  TAUT*  FOUT 

PSI(1,JP,1)  =  PSI(1,JP+F,2) +TAUT* ( U0 ( JP+F , 2 )  +U0 ( JP , 1 ) ) /2 . 

DO  5400  JP  =  NPART-F+1 ,NPART 
J  =  JP  +  Nl *  F 

U0(JP,1)  =  KPOD(l)  -  OMEGA ( 1 ) 

PSI(1,JP,1)  =  PS00(1,J)+  FLOAT ( NPART- JP ) *TAU*U0 ( JP , 1 ) 

DO  5500  Jl=l , NMODE 

A( Jl , 1 )  =  A( Jl , 2 ) +TAUT* ( 55.*F14(J1 )-59 . *F13(J1 )+37. *Fl2(Jl ) 

-  9.*F11( Jl) )/24 . 

THETA ( Jl , 1 )  =  THETA ( Jl , 2 )  +  TAUT* ( 5 5 . * F2 4 ( Jl ) -59 . *F23 ( Jl ) 

+37 . *F22 ( Jl )  -  9 . *F21 ( Jl ) )/24 . 

APP ( Jl , 1 )  =  Fl4 ( Jl ) 

THETA  PRIME  IS  THE  AVERAGE  OF  THE  DISCRETE  AND  FUNCTIONAL 
VALUES  OF  THE  DERIVATIVE 
TP ( Jl , 1 )  =  F24 ( Jl ) 

ATEMP(Jl)  =  19 . *F14 ( Jl )  -  5.*Fl3(Jl)  +  Fl2(Jl) 

ATEMPl ( Jl )  =  19 . *F24 ( Jl )  -  5.*F23(Jl)  +  F22(J1) 

CONTINUE 

DO  5800  M=1 , MAXIT 
DO  5700  Jl=l , NMODE 
CALL  EVOLV (Jl,F25,Fl5,l) 

CORA(Jl)  *  A ( Jl , 1 ) 

A( Jl , 1 )  =  A( Jl , 2 )  +  TAUT* ( 9 . *F15 ( Jl )  +  ATEMP ( Jl ) ) /24 . 

CTHET(Jl)  =  THETA( Jl , 1 ) 

THETA ( Jl , 1 )  =  THETA ( Jl , 2 )  +  TAUT* ( 9 . *F25 ( Jl )  +  ATEMPl ( Jl )) /24 . 
TPC(Jl)  =  F25 ( Jl ) 

TP ( Jl , 1  )  =  F2  5 ( Jl ) 

APP ( Jl , 1 )  -  Fl 5 ( Jl ) 

I F (  MOD ( N , 1 0 )  .EQ.  0  ) THEN 
NPLOT  =  INT( N/10 )  +  1 
TIME ( NPLOT )  =  TIM 
EWAV( NPLOT)  =  0. 

END  IF 

DO  5600  JP  =1 , NPART-4  *F 
CALL  F4 ( JP , FOUT , 1 ) 

U0(JP,1)  =  U0 ( JP+F , 2 )  +TAUT* ( 9 . *FOUT  +TEMP( JP ) )/24 . 

PSI(1,JP,1)  -  PSI ( 1 , JP+F , 2 )  +TAUT* (9.*U0(JP,1) + TEMPI ( JP ) )/24 . 
DO  5750  Jl=l, NMODE 

I F (  ABS (  A( Jl , 1 ) -CORA ( Jl ) ) /ABS ( CORA ( Jl ) )  .GT.  ERROR  .OR. 

ABS ( THETA ( Jl , 1 ) -CTHET ( J 1 ) ) /ABS  ( CTHET ( J 1 ) )  .GT.  ERROR  .OR. 

ABS (TP(Jl,l)— TPC(Jl) ) /ABS ( TP ( J 1 , 1 ) )  . GT .  ERROR2 ) THEN 

GO  TO  5799 
ELSE 

I F (  MOD( N, 10 )  .EQ.  0 ) THEN 
NPLOT  =  I NT ( N/l 0 )  +  1 

PLAR ( NPLOT , J 1 )  =  A(Jl,l)*SIN( THETA ( Jl , 1 ) ) 

PLAI ( NPLOT, Jl )  =  -A( Jl , 1 ) *COS ( THETA ( Jl , 1 ) ) 

FREQ( NPLOT, Jl )  =  TP(J1,1) 

EWAV( NPLOT)  =  EWAV( NPLOT)  +  KWIGR*  *2* ( ( BETAZ 0* OMEGA ( Jl ) * 

A( Jl , 1 ) ) **2  +  KPOD( Jl)**2*(PHIl(Jl,l)**2+PHl2(Jl,l)**2) 

/4 . )/( 4 . *NU* ( GAM0-1 . ) *KWIGL*  *2  ) 

END  IF 


END  IF 

1750  CONTINUE 

GO  TO  5850 

5799  I  F  ( M  .EQ.  MAXIT)WRITE( 6 , 1 )N,M 

>800  CONTINUE 

850  DO  6200  K=4,l,-1 

DO  5900  Jl=l , NMODE 
A(  Jl , K+l )  =  A(Jl,K) 

APP{ Jl , K+l )  =  APP  ( Jl ,  K ) 

THETA( Jl , K+l )  =  THETA ( Jl , K ) 

5900  TP( Jl , K+l )  =  TP ( Jl , K ) 

DO  6100  JP=1 , NPART-F 
U0 ( JP+F , K+l )  =  U0 ( JP , K ) 

PSI(1,JP+F, K+l )  =  PS I ( 1 , JP , K ) 

6100  CONTINUE 

'200  CONTINUE 

DO  6300  Jl=l , NMODE 
Fll(Jl)  =  Fl 2 ( Jl  ) 

F12 ( Jl )  =  F13 ( Jl ) 

Fl 3 ( Jl )  =  Fl 4 ( Jl ) 

F21 ( Jl )  =  F22 { Jl ) 

F22 ( Jl )  =  F23 ( Jl ) 

'300  F23 ( Jl )  -  F24 ( Jl ) 

CALL  SECOND (CPU) 

I F (  CPU  .GE.  . 9 5* FLOAT ( TLIM)  )  GO  TO  7001 
7000  CONTINUE 

301  NLAST  =  N 

JJ  =  0 

DO  7002  JP»5 , NPART ,  5 
JJ  =  JJ  +  1 

U0I( JJ, NCOUNT+1 )  =  U0 ( JP , 2 ) 

/ 002  PSI I ( JJ , NCOUNT+1 )  =  PSI(1,JP,2) 

TPLOT( NCOUNT+1 )  =  TIM 

WRITE ( 4 )  TIME , FREQ , EWAV, PLAR , PLAI , KPOD , OMEGA , KWIGR , 

1  NU,GAM0, BETAW, RISE, FILL, REF, EPS, PHASE, 

1  ERROR, BETAZ0, ERROR2, PSII ,U0I ,TPLOT, TLIM, BETA0 

WRITE ( 5 ) NWIG , NPART , F , NMODE , NPLUS , MAXIT , NSEP , NTIMES 
1  ,NCOUNT, NLAST 

WRITE ( 7 )  PSI ,U0, Fll , Fl 2, Fl 3 , F21 , F22 , F23 , A, APP, THETA, 

1  TP , PHI 1 , PHI  2 

CLOSE ( UNIT=1 ) 

CLOSE ( UNIT=2  ) 

CLOSE( UNIT=3 ) 

CLOSE ( UNIT=4  ) 

CLOSE ( UNIT=5 ) 

CLOSE ( UNIT=7 ) 

END 

SUBROUTINE  F4 ( JP , FOUT , MM ) 

!  REAL  BETAZ , BETAZ0 , GAM , KPOD ( 20 ) , OMEGA ( 20 ) 

REAL  PSI(20,3000,5) ,U0(3000,5) , PH 11(20,5) , PHI 2 (20, 5) ,A(20,5) , 

,  1  APP(  20, 5) ,TP(  20, 5) ,THETA( 20, 5) , FOUT 

INTEGER  JP, MM, NMODE 

'  COMMON/BLKl/BETAZO , GAM0 , BETAW, KPOD , OMEGA 

COMMON/BLK5/PS I , U0 , PHI 1 , PHI  2 , A , THETA , APP , TP 

COMMON/BLK4/ALPHA1 , ALPHA2 , NMODE 

BETAZ  =  BETAZ 0* ( U0 ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  SQRT(  ( 1 .  +  ( GAM0  * BETAW ) *  *  2 ) / ( 1 . -BETAZ  *  BETAZ ) ) 

SUMl  =KPOD ( 1 ) * ( PHI  2 ( 1 , MM ) *COS ( PSI ( 1 , JP , MM ) -THETA ( 1 , MM )  ) 

1  -PH I 1 ( 1 , MM )*SIN(PSI(1,JP, MM ) -THETA ( 1 , MM ) ) ) 

SUM2  =( KPOD( 1 )-BETAZ0* BETAZ *( OMEGA ( 1 )+TP( 1 , MM ) ) ) *A( 1 , MM ) * 

1  SIN( PSI ( 1 , JP , MM ) -THETA ( 1 , MM )  ) -BETAZ0* BETAZ  *APP( 1 , MM ) * 

1  COS ( PSI ( 1 , JP , MM ) -THETA ( 1 , MM ) ) 

DO  100  Jl=2, NMODE 

SUM1=SUM1  +KPOD( Jl ) *( PHI2( Jl , MM ) *COS( PSI ( Jl , JP , MM ) -THETA( Jl , MM ) ) 
1  -PHIl( Jl,MM)*SIN(PSI( Jl, JP, MM) -THETA ( Jl,MM) ) ) 

SUM2=SUM2+ ( KPOD( J 1 ) -BETAZ 0* BETAZ* ( OMEGA( J 1 ) +TP ( Jl , MM ) ) ) *A( Jl ,MM) 


1 

1 


*SIN( PSI ( Jl , JP, MM) -THETA { Jl ,MM) ) -BETAZO*BETAZ*APP ( Jl , MM ) * 

COS (PS I ( Jl, JP, MM) -THETA ( J1,MM) ) 

100  CONTINUE 

FOUT  =  ALPHAl*SUMl*( 1 . -BETAZ*BETAZ )/GAM  +  ALPHA2*SUM2/( GAM* 

1  GAM) 

RETURN 

END 

SUBROUTINE  EVOLV( Jl , K2 , Kl , MM ) 

REAL  BETAZO , BETAZ , GAM , KPOD (20) , OMEGA ( 20) 

REAL  A(20,5) ,PHI1<  20, 5) ,PHI 2(20, 5) , PSI (20, 3000,5) ,U0(  3000,5) 

REAL  K2 ( 20 ) , Kl ( 20  ) , TIM, THETA( 20 , 5 ) ,APP(20,5) ,TP(20,5) 

REAL  NUI , NUR 

INTEGER  Jl , MM , N , F , NPART 

COMMON/BLKl/BETAZ  0 , GAM0 , BETAW, KPOD , OMEGA 
COMMON/BLK5/PSI , U0 , PHI 1 , PHI  2 , A , THETA , APP , TP 
COMMON/BLK 3/TIM , NPART , N , RI SE , TAU , BETAl , BETA2 , F , NUI , NUR 
DUMl  =  0. 

DUM2  =  0. 

DUM3  =  0. 

DUM4  =  0. 

DO  100  JP  =1, NPART 
J  =  JP  +  ( N-l ) *F 

TRISE  =  1.  -  EXP ( -FLOAT ( J-l ) *TAU/RI SE ) 

BETAZ  =  BETAZO  * ( U0 ( JP , MM )  +  OMEGA ( 1 )  )/KPOD(l) 

GAM  =  SQRT(  (l.+(GAM0*BETAW)**2)/(l.-BETAZ*BETAZ)  ) 

PSI  (  Jl , JP,MM)=KPOD( Jl ) * ( PSI ( 1 , JP,MM)+OMEGA( 1 ) *TIM)/KPOD( 1 ) 

1  -  OMEGA ( Jl ) *TIM 

DUMl  =  DUMl  +  COS(PST(Jl,JP, MM ) -THETA ( Jl , MM ) ) *TAU*TRI SE*GAM0/GAM 
DUM2  =  DUM2  +  S IN ( PS I ( Jl , JP , MM ) -THETA( Jl , MM ) ) *TAU*TRI SE*GAM0/GAM 
DUM3  =  DUM3  +  COS ( PS I ( Jl , JP , MM ) -THETA( Jl , MM ) ) *TAU*TRI SE 
100  DUM4  =  DUM4  +  SIN ( PSI ( Jl , JP , MM ) -THETA( Jl , MM ) ) *TAU*TRI SE 

K1(J1)  =  -NUR*A( Jl , MM )/2 . -  BETAl *DUM2/OMEGA( Jl ) 

K2(J1)  =  -NUI/2 .  +  BETAl *DUMl/( A( J 1 , MM ) * OMEGA ( Jl ) ) 

PHI 1 ( Jl , MM )  *  -  BETA2*DUM3/KPOD( Jl  )**2 
PHI 2 ( Jl , MM )  *  -  BETA2 *DUM4/KPOD ( Jl ) * *2 
RETURN 
END 


JJJJJJJJJJ 

iJJJJJJJJJ 
JJJJJJJJJ 


9999999999999999999999999999999999999999999999999999 
Digital  Equipment  Corporation  -  VAX/VMS  Version  V4.6 
9999999999999999999999999999999999999999999999999999 


JJJJJJJJJJ 

JJJJJJJJJJ 

JJJJJJJJJJ 


M 

M 

AAA 

RRRR 

AAA 

BBBB 

L 

EEEEE 

MM 

MM 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M  M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

RRRR 

A 

A 

BBBB 

L 

EEEE 

M 

M 

AAAAA 

R  R 

AAAAA 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

B 

B 

L 

E 

M 

M 

A 

A 

R 

R 

A 

A 

BBBB 

LLLLL 

EEEEE 

SSSSSSSS 

SSSSSSSS 


PPPPPPPP 

PPPPPPPP 


RRRRRRRR 

RRRRRRRR 


AAAAAA 

AAAAAA 


DDDDDDDD 

DDDDDDDD 


BBBBBBBB 

BBBBBBBB 


SS 

pp  PP 

RR 

RR 

AA 

AA 

DD 

DD 

BB  BB 

ss 

PP  PP 

RR 

RR 

AA 

AA 

DD 

DD 

BB  BB 

SS 

pp  PP 

RR 

RR 

AA 

AA 

DD 

DD 

BB  BB 

ss 

pp  pp 

RR 

RR 

AA 

AA 

DD 

DD 

BB  BB 

ssssss 

PPPPPPPP 

RRRRRRRR 

AA 

AA 

DD 

DD 

BBBBBBBB 

ssssss 

PPPPPPPP 

RRRRRRRR 

AA 

AA 

DD 

DD 

BBBBBBBB 

SS 

pp 

RR 

RR 

AAAAAAAAAA 

DD 

DD 

BB  BB 

SS 

pp 

RR 

RR 

AAAAAAAAAA 

DD 

DD 

BB  BB 

ss 

pp 

RR 

RR 

AA 

AA 

DD 

DD 

BB  BB 

ss 

pp 

RR 

RR 

AA 

AA 

DD 

DD 

BB  BB 

SSSSSSSS 

pp 

RR 

RR 

AA 

AA 

DDDDDDDI 

) 

BBBBBBBB 

SSSSSSSS 

pp 

RR 

RR 

AA 

AA 

DDDDDDDI 

) 

BBBBBBBB 

FFFFFFFFFF 

OOOOOO 

RRRRRRRR 

•  • 

9  9 

11 

FFFFFFFFFF 

000000 

RRRRRRRR 

•  • 
t  9 

11 

FF 

00 

00 

RR 

RR 

9  9 

1111 

FF 

00 

00 

RR 

RR 

•  • 

9  9 

1111 

FF 

00 

OO 

RR 

RR 

11 

FF 

00 

00 

RR 

RR 

11 

FFFFFFFF 

00 

00 

RRRRRRRR 

•  • 

9  9 

11 

FFFFFFFF 

00 

00 

RRRRRRRR 

9  9 

11 

FF 

00 

00 

RR 

RR 

•  • 

9  9 

11 

FF 

00 

00 

RR 

RR 

•  • 

9  9 

11 

,  ,  ,  , 

FF 

00 

00 

RR 

RR 

11 

•  •  •  • 

FF 

00 

00 

RR 

RR 

11 

•  •  •  • 

FF 

000000 

RR 

RR 

9  9 

min 

,  ,  ,  , 

FF 

000000 

RR 

RR 

9  9 

min 

File  _VC$DRBl : [ MARABLE . TWSTAG J SPRADB . FOR ; 1  (4143,1,0),  last  revised  on 

14-DEC-1984  10:33,  is  a  10  block  sequential  file  owned  by  UIC  [MARABLE].  The 
ecords  are  variable  length  with  a  fixed  control  size  of  2  bytes  and  implied 
(CR)  carriage  control.  The  longest  record  is  71  bytes. 

lob  SPRADB  (1985)  queued  to  LN03QUE  on  21-MAR-1988  13:48  by  user  MARABLE,  UIC 
[MARABLE],  under  account  4790  at  priority  100,  started  on  printer  LTA7 :  on 
21-MAR-1988  13:48  from  queue  VCLN03A. 
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TEST  PROGRAM  FOR  MODEL  EQUATIONS  IN  HIGH  POWER 
FEL  PROBLEM. 

REVISION  12/6  NEW  FORMULATION  OF  WAVE  EQUATION 
REVISION  12/14  ADDITION  OF  DOUBLE  PRECISION 
REAL  A(  1200) , EPS, TAU, OMEGA 
REAL  TIME( 1200 ), KPOD,GAM( 400 ) 

REAL  GROWTH ( 1200 ) , KWIGR , KWIGL , NU, FEQ( 1200 ) ,RISE 

DOUBLE  PRECISION  PSIM (  4 00 ) , PS I  0 ( 4 00 ) , BETAl , DUM , DUMl , PS0 0 

DOUBLE  PRECISION  UOLD ( 400 ) , UNEW ( 400 ) , FEQD , AD , GROW , THETAP 

DOUBLE  PRECISION  ALPHAl ,AMAG , THETA 

INTEGER  J , NPART, NTIMES , N, F 

CHARACTER* 40  GLAB , XLAB , YLAB 

CHARACTER* 80  ABl , AB2 , AB3 

PS00(J)  =  -OMEGA*DFLOAT( J)*TAU 

OPEN ( UN IT* 2 , NAME= ' OUT . DAT ' , STATUS* ' NEW' ) 

GLAB*  'WAVE  AMPLITUDE  VS.  TIME$ ' 

XLAB  ='  TIME$ ' 

YLAB  =  '  WAVE  AMPLITUDE$ ' 

WRITE( 6,101) 

FORMAT ( '  INPUT  NO.  OF  PART., NO.  OF  INTERAT . , GAM , BETAWIG , 

1  KWIGr ,BUDKFR, KPOD, OMEGA, KWIGL, EPS, F, RISE  ') 

READ ( 5 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , KPOD , OMEGA , KWIGL , EPS  ,  F 

WRITE( 6,103) 

WRITE( 2,103) 

FORMAT ( '  THE  INPUT  DATA: NPART , NTIMES , GAM , BETAW, KWIGR , NU , KPOD 
1  , OMEGA, KWIGL, EPS, F, RISE  IS:  ') 

WRITE ( 6 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , KPOD , OMEGA , KWIGL , EPS  , 

WRI TE ( 2 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , KPOD , OMEGA , KWIGL , EPS  , 

FORMAT ( ' NPART* ' ,I4,3X, ' NTIMES= ' ,I4,3X, ' GAM* ' ,F8.4,3X, 

1  ' BETAWIG= ' , F8 . 4 ) 

FORMAT (  ' KWIGR* ' ,F8.4,3X, ' NU= ' ,E10.4,3X,  'KPOD*' , El  0.4 
1 , 3X , ' OMEGA* ' , El  0.4) 

FORMAT ( ' KWIGL* ' , F8 .4 , 3X, 'EPS*'  , E10 . 4 , 3X , ' F= '  , 1 3 , 3X , 

1 ' RISE* ' ,F8.4) 

WRITE ( ABl ,104) NPART , NTIMES , GAM0 , BETAW 
WRTTF.r  AR2 , 105  )  KWIGR,  NU  ,  KPOD,  OMEGA 
WRITE(AB3, 106 (KWIGL, EPS, F, RISE 
BETA0  =  SQRT(1.  -1 ./( GAM 0* GAM 0 ) ) 

BETAZ0  =  SQRT( BETA0 *BETA0  -  BETAW*BETAW/2 .  ) 

BETAl  =  4 . *NU* BETAW* ( KWIGL ) *  *  2/( BETAZ 0*BETA0* GAM 0  * OMEGA 
1  * ( KWIGR*KWIGR ) ) 

ALPHAl  =  -BETAW/( 2 . *BETAZO*BETAZ0 ) 

INITIALIZE  PHASE  AND  AMPLITUDE 
AMAG  =  EPS 
THETA  =  0. 

A( 1 )  -  EPS 
TIME ( 1  )  =  1. 

TAU  =1 ./FLOAT ( NPART- 1 ) 

DO  100  J=l, NPART 

PSl0(J)  =  PS00(J-1)  +  (KPOD  -  OMEGA) * (NPART  -  J)*TAU 
UOLD(J)  =  KPOD  -  OMEGA 
CONTINUE 
WRITE( 6,102) 

WRITE ( 2 , 102 ) 

FORMAT ( '  THE  WAVE  AMPLITUDES  ARE:  ') 

BEGIN  LOOP  FOR  TIME  INCREMENTS 
DO  1000  N  =2, NTIMES 

TIME ( N )  =  1.  +  FLOAT ( N-l ) * FLOAT ( F )/( NPART  -1) 

BEGIN  PART  II:  STEP  AMPLITUDES  AND  PHASES 
COMPLETE  SUM  FOR  AMPLITUDE  STEP 
DUM  =  0. 

DUMl  =  0. 

DO  200  J=l, NPART 
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BETAZ  =  BETAZO* ( UOLD( J )  +  OMEGA ) /KPOD 

GAM ( J ) =  SQRT(  (l.-BETAZ*BETAZ)/(l./(GAM0*GAM0)+BETAW*BETAW/2. )  ) 

DUM  =  DUM  +  DCOS ( PSI 0 ( J )  -  THETA ) *TAU*GAM ( J ) 

DUMl  *  DUM1  +  DS IN( PSI 0 { J )  -  THETA ) *TAU*GAM ( J ) 

CONTINUE 

GROW  =  BETAl * ( 1 . -EXP ( ( 1 . -TIME ( N ) ) /RISE ) ) * DUM/ AMAG 
THETAP  =  BETAl * ( 1 . —EXP ( ( 1 . -TIME ( N ) ) /RISE ) ) * DUM 1 /AMAG 
AD  =  AMAG  +  D FLOAT ( F ) *TAU* GROW* AMAG 
A(N)  =  AD 
DO  300  J  =  1  , NPART 

PSIM(J)  =  PSIO(J)  +  D FLOAT ( F ) *TAU*UOLD ( J ) 

UNEW(J)  =  UOLD(J)  +  DFLOAT ( F ) *TAU*ALPHAl *GAM ( J ) *GAM ( J ) * ( 

1  ( KPOD* KPOD-BETAZO* BETAZ 0  *OMEGA* ( UOLD ( J ) +OMEGA )  )* 

1  ( AMAG* DCOS ( PSI 0 ( J )  -  THETA))  + 

1  BETAZO* BETAZO* (UOLD( J )+OMEGA) * ( AMAG*GROW*DSIN( PSIO ( J ) -THETA) 

1  -  AMAG* THETAP* DCOS ( PS 1 0 ( J ) -THETA )  )) 

CONTINUE 

END  OF  PART  II  A) 

BEGIN  BOOKEEPING 
DO  400  J=1 , NPART-F 
PSIO(J)  =  PSIM( J  +  F ) 

UOLD(J)  =  UNEW( J+F ) 

CONTINUE 
DO  500  J=F , 1 , -1 

UOLD ( NPART+1— J ) =KPOD-OMEGA+ ( J-l ) *TAU*ALPHAl*  ( 

1  ( KPOD*KPOD-BETAZO*BETAZO*OMEGA*KPOD) * 

1 ( AMAG*DCOS ( PS00 ( NPART- J+ ( N-l ) *  F ) -THETA )  ) 

1 +BETAZ0*BETAZ0* KPOD* (AMAG*GROW*DS IN (PS00( NPART- J+(N-1 ) *F) -THETA 

1  -  AMAG*THETAP*DCOS ( PS00 ( NPART-J+ ( N-l ) *F ) -THETA )  )) 

PSIO ( NPART+1— J ) =PS00 ( NPART- J+ ( N-l ) *F ) + ( J-l ) *TAU* ( KPOD-OMEGA 
1  +UOLD( NPART+1— J )  )/2. 

CONTINUE 

END  OF  PART  II  B)  BOOKEEPING 
AMAG  =  AD 

GROWTH (N-l)  =  GROW 

THETA  =  THETA  +  DFLOAT( F ) *TAU*THETAP 
FEQD  =  OMEGA  +  THETAP 
FEQ(N-l)  =  FEQD 

REPEAT  PHASE  AND  AMPLITUDE  INCREMENT 
FOR  NEXT  TIME  STEP 

I F (  FLOAT ( N ) /I 0 .  -  N/10  .EQ.  0.)  THEN 
WRITE (  6  ,  *  )  A ( N ) , TIME ( N ) 

WRITE (  2  ,  *  )  A ( N ) , TIME ( N ) 

END  IF 
CONTINUE 

GROWTH ( NTIMES )  =  GROW 
FEQ(NTIMES)  =  OMEGA  +  THETAP 
WRITE ( 6 ,1001 ) 

WRITE ( 2 ,1001 ) 

FORMAT ( '  THE  FINAL  PHASES  ARE:') 

WRITE (6,*) ( PSIM( J ) ,J=1, NPART) 

WRITE ( 2 , * ) ( PSIM( J) , J=1 ,NP^RT) 

CALL  ANOTAT ( %REF ( XLAB ) , %REF ( YLAB ) ,1,0,0,DSHL) 

CALL  PWR I T ( 500,80, %REF( ABl ) ,70, 1 ,0,0 ) 

CALL  PWRIT( 500,48, %REF( AB2 ) ,72,1,0,0) 

CALL  PWRIT( 500,16, %REF( AB3 ) ,75,1,0,0) 

CALL  EZXY(TIME, A, NTIMES, %REF(GLAB) ) 

YLAB  =  '  GROWTH  RATE  $ ' 

GLAB  =  '  GROWTH  RATE  NORMALIZED  TO  TRANSIT  TIME$ ' 

CALL  ANOTAT ( %REF( XLAB) , %REF ( YLAB ) , 1,0,0, DSHL) 

CALL  EZXY(TIME, GROWTH, NTIMES, %REF(GLAB)  ) 

YLAB  =  '  REAL  FREQUENCY  $ ' 

GLAB  =  '  REAL  FREQUENCY  VS.  TIME  $' 

CALL  ANOTAT ( %REF ( XLAB ) , %REF ( YLAB ) ,1,0,0, DSHL) 

CALL  EZXY ( TIME, FEQ, NTIMES, %REF( GLAB) ) 
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TEST  PROGRAM  FOR  MODEL  EQUATIONS  IN  HIGH  POWER 
FEL  PROBLEM. 

REVISION  12/3  TO  CORRECT  BOUNDARY  CONDITIONS 

REVISION  1/16  TO  INCLUDE  PHENOMENOLOGICAL  DAMPING 

REAL  A( 3000 ) ,PSIM(4000) ,PSI0(4000) , EPS , BETAl , TAU , DUM , OMEGA 

REAL  TIME( 3000 ) , KrOD, UOLD( 4000 ) ,UNEW(4000) ,GAM(4000) ,DUMl 

REAL  GROWTH ( 3000 ) , KWIGR , KWIGL , NU, RISE , GROW2 (3000 ) 

INTEGER  J , NPART , NTI MES ,  N ,  F 
PARAMETER  ( PI=3 . 141592653589 ) 

CHARACTER* 40  GLAB , XLAB , YLAB 
CHARACTER* 80  ABl , AB2 , AB3 , AB4 
PS00(J)  =  -OMEGA* J*TAU 

OPEN ( UNIT=2 , NAME* ' OUT . DAT ' , STATUS* ' NEW ' ) 

FORMAT ('WAVE  AMPLITUDE  W/  REFL=',F5.3) 

XLAB  =  '  TIME$ ' 

YLAB  =  '  WAVE  AMPLITUDE$ ' 

WRITE( 6,101) 

FORMAT ( '  INPUT  NO.  OF  PART., NO.  OF  INTERAT . , GAM , BETAWIG , 

1  KWIGr , BUDKER , NWIG, EPS, F, RISE, NPLUS , REF  ') 

READ ( 5 , * ) NPART , NTIMES , GAM0 , BETAW, KWIGR , NU , NWIG , EPS , F , RI SE , NPLUS , 

WRITE (6,103) 

WRITE ( 2 , 103 ) 

FORMAT ( '  INPUT  DATA: NPART , NTIMES , GAM , BETAW , KWIGR , NU , NWIG 
1, EPS, F, RISE, NPLUS, REF  IS:') 

WRITE ( 6 , * ) NPART , NTIMES , GAM0 , BETAW , KWIGR , NU , NWIG , EPS , F , RISE , NPLUS 

WRITE ( 2, * ) NPART, NTIMES, GAM0 , BETAW, KWIGR, NU, NWIG, EPS, F, RISE, NPLUS 

KWIGL  =  2 . *FLOAT(NWIG) *PI 

BETA0  =  SQRT(1.  -  1 ./( GAM0*GAM0 ) ) 

BETAZ0  -  SQRT( BETA0*BETA0  -  BETAW*BETAW/2 . ) 

NOPT  -  JINT( 2 . * FLOAT ( NWIG) *BETAZ0/( 1 . — BETAZ0 ) )  +NPLUS 
OMEGA  =  ( FLOAT ( NOPT ) *PI ) /BETAZ0 
KPOD  =  KWIGL  +  BETAZ  0  *OMEGA 
WRITE ( AB4 , 1 ) REF 

WRI TE ( ABl , 1 0  4 ) NPART , NTIMES , GAM0 , BETAW , NOPT 

FORMAT ( ' NPART* ' ,I4,2X, ' NTIMES= ' ,I4,2X, ' GAM= ' ,F8.4,2X, 

1  ' BETAWIG= ' , F8 . 4 , 2X , ' NOPT= ',14) 

FORMAT ( ' KWIGR= ' , F8.4.3X, 'NU=' ,E10.4,3X, 'KPOD*' , El  0.4 
1,3X, 'OMEGA*' , El 0.4) 

WRITE ( AB2 , 1 0  5 ) KWIGR , NU , KPOD , OMEGA 

WRITE(AB3,106)KWIGL,EPS,F, RISE, NWIG, NPLUS 

FORMAT ( 'KWIGL*'  ,F8.4,2X, 'EPS*' , El  0 . 4 , 2X , ' F= ' , I  3 , 2X , 

1 ' RISE* ' ,F8.4,2X, 'NWIG*' ,I3,1X, 'NP=' ,13) 

BETAl  =  4 . *NU* BETAW* ( KWIGL ) *  *  2/( BETAZ 0  *BETA0*GAM0  *OMEGA 
1  * ( KWIGR*KWIGR ) ) 

ALPHAl  =  -BETAW/( 2 . *BETAZ0*BETAZ0 ) 

ALPHAR2  =  (1.  -  REF ) /BETAZ0 

ALPHAl 2  =-4 . *NU*KWIGL**2* ( 1 . -BETAW* *  2/2 . )/( BETAZ0**2*KWIGR**2 
1  * GAM 0* OMEGA) 

INITIALIZE  PHASE  AND  AMPLITUDE 
ARPRIM  =  0. 

AIMPRIM  =  0. 

A ( 1 )  =  EPS 
TIME ( 1 )  =  1. 

TAU  =1 ./FLOAT (NPART-1 ) 

TAUT  =  FLOAT ( F ) *TAU 
ARN  =  EPS*COS ( 3 . *TAUT ) 

AIMN  *  EPS*SIN( 3 . *TAUT ) 

DO  100  J=l, NPART 

PSI0(J)  *  PS00(J-1)  +  (KPOD  -  OMEGA ) * ( NPART  -  J)*TAU 
UOLD(J)  =  KPOD  -  OMEGA 
CONTINUE 
WRITE ( 6 , 102 ) 

WRITE ( 2 , 102 ) 
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FORMAT ( '  THE  WAVE  AMPLITUDES  ARE:  ' ) 

BEGIN  LOOP  FOR  TIME  INCREMENTS 
GROWTH ( 1 )  =  0. 

GROW 2 ( 1 )  =  0 . 

DO  1000  N  =2 , NTIMES 

TIME ( N )  =  1.  +  FLOAT ( N-l ) *TAUT 

TRISE  =  (1.  -EXP (  ( 1 .-TIME(N) )/RISE  )) 

BEGIN  PART  II:  STEP  AMPLITUDES  AND  PHASES 
COMPLETE  SUM  FOR  AMPLITUDE  STEP 
DUM  =  0. 

DUMl  =  0. 

DO  200  J  =  1 , NPART 

BETAZ  =  BETAZ0* ( UOLD( J )  +  OMEGA ) /KPOD 

GAM ( J ) =  SQRT(  (l.-BETAZ*BETAZ)/(l./(GAM0*GAM0)+BETAW*BETAW/2. )  ) 

DUM  =  DUM  +  COS ( PSI 0 ( J ) ) *TAU*GAM( J ) 

DUMl  =  DUMl  +  SIN(PSI0( J) )*TAU*GAM( J) 

CONTINUE 

ARNl  =  ARN  +TAUT* ( ( ALPHAI2*AIMN+BETA1*DUM ) *TRI SE-ALPHAR2 *ARN ) 
AIMNl=AIMN  +TAUT* { ( -ALP HA 1 2  *ARN+BETAl *DUMl ) *TRI SE-ALPHAR2  *AIMN ) 
A(N)  =  SQRT( ARNl* ARNl  +  AIMNl *AIMNl ) 

DO  300  J=l, NPART 

PSIM(J)  =  PSI0(J)  +  TAUT*UOLD( J ) 

UNEW(J)  =  UOLD(J)  +  TAUT*ALPHAl *GAM ( J ) *GAM( J ) * ( 

1  ( KPOD*KPOD-BETAZO*BETAZO*OMEGA* ( UOLD( J ) +OMEGA) )* 

1  ( ARN* COS ( PSI 0 ( J ) ) +AIMN*SIN( PSI 0  (  J )  ) )  + 

1  BETAZ 0* BETAZ 0* ( UOLD( J ) +OMEGA ) * ( ARPRIM*SIN( PSI0 ( J ) ) 

1  -  AIMPRIM*COS ( PSI0 ( J ) )  )) 

CONTINUE 

END  OF  PART  II  A) 

BEGIN  BOOKEEPING 
DO  400  J=1 , NPART— F 
PSI0(J)  =  PSIM( J+F ) 

UOLD(J)  =  UNEW( J+F ) 

CONTINUE 
DO  500  J=F ,1,-1 

UOLD ( NPART+1- J ) -KPOD-OMEGA+ ( J-l ) *TAU*ALPHAl * ( 

1  ( KPOD*KPOD-BETAZ0*BETAZ0*OMEGA*KPOD) * 

1 ( ARN* COS ( PS00 ( NPART- J+ ( N-l ) *F ) ) +AIMN* SIN ( PS00 ( NPART- J +( N-l )* F ) ) 

1  +BETAZ 0* BETAZ 0* KPOD* ( ARPRIM*SIN ( PS00 ( NPART-J+ (N-l ) *F ) ) 

1  -  AIMPRIM*COS ( PS00 ( NPART- J+ ( N-l ) * F ) ) )  ) 

PS I 0 ( NPART+1- J ) =PS00 ( NPART- J+ (N-l ) *F)+( J-l ) *TAU* ( KPOD-OMEGA 
1  +UOLD( NPART+l-J )  )/2 . 

CONTINUE 

END  OF  PART  II  B)  BOOKEEPING 

ARPRIM  -  ( ALPHAI 2  *AIMN+BETAl  *DUM  )  *TRI  SE-ALPHAR2  *ARN 
AIMPRIM  =  ( -ALPHAI 2  *ARN+BETAl *DUMl ) *TRI SE-ALPHAR2  *AIMN 
GROWTH ( N )  = ( ARNl * ARPRIM  +  AIMNl *AIMPRIM )/( A( N ) *A( N ) ) 

GROW2 ( N )  =  (A(N)  -  A( N-l ) )/( A( N-l ) *TAUT ) 

ARN  =  ARNl 
AIMN  =  AIMNl 

REPEAT  PHASE  AND  AMPLITUDE  INCREMENT 
FOR  NEXT  TIME  STEP 

I F (  FLOAT ( N ) /I 0 .  -  N/10  .EQ.  0.)  THEN 
WRITE ( 6 , * )  A ( N ) , TIME ( N ) 

WRITE ( 2 , * )  A( N ) , TIME ( N ) 

END  IF 
CONTINUE 
WRITE( 6,1001 ) 

WRITE ( 2 , 1001 ) 

FORMAT ( '  THE  FINAL  PHASES  ARE:') 

WRITE (6,*) ( PSIM( J ) , J-l, NPART) 

WRITE ( 2, * ) ( PSIM( J) , J=1 , NPART) 

CALL  ANOTAT ( %  REF ( XLAB ) , %REF ( YLAB ) ,1,0,0, DSHL ) 

CALL  PWRIT( 500,80, %REF ( ABl ) , 70 , 1 , 0 , 0 ) 

CALL  PWRIT(500,48,%REF(AB2) ,72,1,0,0) 


CALL  PWRIT( 500,16,  %REF( AB3 ) ,75,1,0,0) 

CALL  EZXY ( TIME , A , NTIMES , %REF ( AB4 ) ) 

YLAB  =  '  GROWTH  RATE  $' 

GLAB  =  '  GROWTH  RATE  NORMALIZED  TO  TRANSIT  TIME$ ' 
CALL  ANOTAT ( %REF ( XLAB ) , %REF( YLAB) ,1,0, 0 , DSHL ) 

CALL  EZXY(TIME, GROWTH, NTIMES, %REF(GLAB)  ) 

GLAB  =  'GROWTH  RATE  EVALUATED  BY  DERIV.$' 

CALL  ANOTAT ( % REF (XLAB) , %REF ( YLAB ) ,1,0,0, DSHL) 

CALL  EZXY (TIME, GROW2 , NTIMFS , %REF ( GLAB ) ) 

STOP 

END 
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HIGH  GAIN  FREE  ELECTRON  LASER  OSCILLATORS 


I.  Introduction 

We  have  conducted  an  analytical  and  numerical  analysis  of  the  field  evolution  in  a  high 
gain  free  electron  laser  operating  in  the  oscillator  configuration,  as  depicted  in  Fig.  1.  The 
analysis  is  applicable  to  systems  with  electron  beam  pulse  lengths  which  are  longer  that  the 
particle  transit  time  in  the  resonator.  The  electron  beam  equilibrium  is  therefore  assumed 
to  be  spatially  uniform  and  temporally  stationary.  The  radiation  field  and  phase  averages 
which  are  performed  with  the  ensemble  of  electrons  is  conducted  for  an  interaction  length 
which  consists  of  the  entire  wiggler  structure.  This  is  in  contrast  to  other  simulations 
(theories)  which  perform  the  ensemble  average  over  the  wavelength  of  the  ponderomotive 
potential;  as  is  applicable  to  systems  with  temporally  stationary  fields1  or  short  beam 
pulses2  that  are  spatially  periodic. 

We  find  that  the  numerical  simulations  yield  qualitative  and  quantitative  agreement 
with  the  theory.  The  theory  for  the  example  given  (strong  pump  Compton  regime)  can 
be  separated  into  three  operation  regimes  which  we  shall  denote  as  the  ultra-high  gain, 
moderate  gain  and  low  gain  regimes.  Both  the  ultra-high  gain  {T^L  >>  1)  and  the  low 
gain  (r iqL  «  1)  regimes  yield  growth  rates  that  exhibit  the  same  scaling  with  beam 
current,  energy  and  wiggler  field  as  is  obtained  for  an  FEL  amplifier  operating  in  these 
regimes.  Additionally,  we  consider  a  moderate  gain  regime  ( T^L  >  1)  which  is  of  direct 
interest  to  NRL  experimental  parameters. 

II.  Theoretical  Model 

An  analysis  of  the  space  time  evolution  of  the  fields  and  particles  within  an  FEL  os¬ 
cillator  requires  a  self-consistent  coupling  of  the  fundamental  equations  for  the  particles 
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and  fields^5-^.  We  have  considered  a  Maxwell-Vlasov  description  of  the  fields  and  par¬ 
ticles.  The  analysis  in  Appendix  A  results  in  the  following  system  of  equations  for  the 
coupling  of  the  fields  and  particles.  The  backward  travelling  wave  evolves  according  to, 
dd(,{z)/dz  —  iaaifz)  =  0.  The  forward  travelling  potential  and  the  electrostatic  potential 
evolve  according  to, 

Q 

—af{z)  +  iaaf{z)  = 

-Ci  f  dz,(z'  -  z)expf-iAA'(r'  -  z)](d/3z'  -  iK)[0wdf{z')/?  -  <p{z')} 

Jo 

+  ic3  [*  dz'exp[-iAK{z'  -  z)](d/dz'  -  »JT)|(1  +  fi0)Mf{z,)J 2  -  *(*')],  (l) 

Jo 

j;6(z)  -  iK${z)/2  = 

-d  [  dz’{z'  -  z)exp\-iAK{z'  -  z)\(d/dz'  -  iK)\0v,dj{z')l2  -  <£(z')J 
Jo 

+  ic<  P  dz'  exp[-iAK{z’  -  z)\(d/dz'  -  iK){0waf(z')/2  -  (l  -  0]o)o{z')\. ,  (2) 
Jo 

The  parameters  in  Eqs.  (1)  and  (2)  are  given  by,  a  =  Au/c-u/£(l-/92  )/2w'oC7,  cx  = 
c2(l  -  /5?0)(cj0  +  Au)/vz0,  c2  -  ojp 0m /2ui0 c$ZQ y ,  cs  =  CiUofd^Kc,  c4  =  c3/{  1  -  0%), 
K  =  k0  +  km  and  Aif  —  K  -  (w0  +  Aw)/vjo.  By  making  use  of  the  convolution  theorem, 
the  Laplace  transform  of  Eqs.  (l)  and  (2)  yields, 

{(a  +  *ot)[(3  -  iAJT)2+2cj]  +  iKc1,3w/2}df(s)  = 

{fa  -  iAKf  +  2cs(l  -  0l/A)}af(O),  (3) 

{(a  -  j'AiC)2  +  2cj  }o(a)  =  Cj^fd/fs)  -  iaf(Q)/K\.  (4) 

where  d/(a)  and  o(a)  are  the  Laplace  transformed  vector  and  scalar  potentials,  and  we  have 
retained  only  the  terms  in  the  driving  current  that  arise  from  the  momentum  derivatives  of 
the  phase,  3  A Kjdpz-  We  have  also  assumed  that  the  electron  beam  enters  the  resonator 
unbunched  so  that  <p(z  =  0)  =  0. 

Since  the  singularities  of  the  Laplace  transformed  potentials  are  isolated  poles,  the 
Bromwich  inversion  of  these  transforms  can  be  easily  performed. 
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(5) 


$(2)  =  52Residue{^(i)’ai}exP(JJ:r)> 

3 

af(z)  =  J^^siduefa/fa),  jy}exp(s^),  (6) 

3 

where  ,  are  the  poles  of  the  Laplace  transformed  potentials.  The  solution  for  the  radiation 
potential  given  in  Eq.  (6)  is  of  the  same  generic  form  as  the  solution  obtained  by  Bernstein 
and  Hirshfield9  for  the  FEL  amplifier  configuration.  Our  analysis  shall  differ  in  that  the 
backward  travelling  waves  are  not  neglected,  and  the  combination  of  the  forward  and 
backward  travelling  waves  are  required  to  satisfy  the  appropriate  boundary  conditions 
at  the  mirror  surfaces.  Specifically,  the  tangential  component  of  the  electric  field  must 
be  zero  at  the  non-transmitting  mirror  surface,  i.e.,  5/(0)  +  5b  (0)  =  0.  At  the  partially 
transmitting  mirror  surface  at  the  end  of  the  resonator,  the  tangential  components  of 
the  electric  field  must  be  continuous,  i.e.,  5/(X)exp(— ik0L)  +  5b(£)  exp(ik0L)  —  (1  — 
y/R)aif(L)  exp(-ikoL)  =  0,  where  L  is  the  length  of  the  resonator  and  R  is  the  fractional 
power  reflected  from  the  far  end  resonator  mirror.  This  yields  the  following  expression  for 
the  boundary  conditions  at  the  mirror  surfaces, 

exp  =  eXp{i(^^(l  -  0H2)  -  2 *oj!}  Residu e{5/(s),  *,•}  exp(s/i),  (7) 

which  is  the  equation  that  seif-consistentiy  determines  the  complex  operating  frequency, 
Aw,  of  the  oscillator. 

HI.  Results  for  Compton  Regime 

In  the  Compton  regime  the  effect  of  the  electrostatic  potential  can  be  neglected^10,11). 
In  addition  we  shall  assume  that  the  spatial  derivative  of  the  vector  potential  in  the  driving 
current  is  negligible.  These  derivatives  are  negligible  when,  J(s  -  j'Aif)2]  >>  2cs  and 
|a|  <<  K.  Under  these  conditions  the  Laplace  transform  of  the  vector  potential  is  given 


(8) 


af(a)  =  5/(0)(«  -  iAK)2/  J]  [a  -  «y), 

j=i 

where  ay  are  the  roots  of  the  dispersion  relation,  (a  -  iAJT)2(a  +  «a)  =  —icl8ulKf2  « 
-iw2(l  +  &zo)&tkv  j4/}*Qc2i.  Since  fc0  is  a  free  parameter,  choose  k0  such  that,  AK  = 

k0  +  kw  -  (w0  +  Aw)/vz0  =  -Aw/c  + w2(l  -  £2 /2)/2w0c7.  Which  results  in  the  following 
solutions  to  the  dispersion  relation, 


_  .rAu  u 
*J  c  2wqC7 


r  »2/v^ 


f=(l-^/2))+r0<  1-  i/y/z 

[  -1  -  iJy/Z 


(9) 


where  T0  =  +  0zo )0l/4klc?7i]l/i /Vzo  is  a  spatial  growth  rate  corresponding 

to  the  largest  spatial  growth  rate  in  the  amplifier  case12.  By  evaluating  the  residues  of 
a/(a)  for  each  of  these  poles,  one  obtains  the  following  solution  for  the  spatial  structure  of 
the  radiation  potential, 


«/(s)  =  H  «P(*y4  U°) 

y=l 

It  is  evident  from  Eq.(lO)  that  the  spatial  growth  of  the  radiation  field  can  be  de¬ 
scribed  by  the  interference  of  three  modes;  which  can  be  identified  as  the  positive  and 
negative  energy  beam  modes,  and  a  transverse  electromagnetic  mode.  The  constructive  or 
destructive  nature  of  this  interference  is  dependent  on  the  values  of  the  physical  parame¬ 
ters  which  characterize  the  roots,  jy.  For  physical  parameters  such  that,  r0T  >>  1,  the 
unstable  mode  dominates  and  one  obtains  exponential  spatial  growth  at  the  rate  IV 
The  temporal  growth  rate  of  the  radiation  field  is  obtained  from  the  negative  imaginary 
part  of  the  complex  oscillator  frequency,  Aw.  The  oscillator  frequency  is  determined  by 
the  boundary  conditions  as  expressed  in  Eq.  (7),  which  for  the  approximate  roots  under 
consideration  yields, 

expif^£  -  ^exp{ij-^(l  -  Jl/2)  -  2k0\L)  ]T  exp(s;X).  (11) 

c  3  2w0C7 
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We  shall  consider  three  distinct  solutions  to  this  equation  for  the  complex  oscillator  fre¬ 
quency.  The  first  of  which  is  the  ultra-high  gain  regime  (ToX  >>  1),  in  which  case,  only 
the  fastest  growing  mode  in  Eq.  (11)  is  retained.  The  second  case  is  the  moderate  gain 
regime  (r0X  >  1),  where  only  the  decaying  mode  in  Eq.  (10)  is  neglected.  The  final  case 
is  valid  for  arbitrary  gain  and  all  terms  in  Eq.  (ll)  are  retained.  The  imaginary  part  of 
the  oscillator  frequency  yields  the  following  temporal  growth  rates, 

rw-  =  ^ In  +  ToX/2,  Ultra-High  Gain  (12) 

C  m  3 

ru-  =  iln^  +  j  ln[2  cosh(r0X)  +  2  cos(V§r0Z)], 

C  M  V 

Moderate  Gain  (13) 

Tu-  =  ^ln^^  +  jln[l  +  4cos(\/3roI)cosh(roZ)  +  4cosh2(ro£)]. 

Arbitrary  Gain  (14) 

In  each  of  the  expressions  for  the  growth  rate  the  first  term  is  negative  definite.  This 
represents  the  effect  of  losses  at  the  mirrors  and  the  coupling  losses  due  to  the  splitting  of 
the  radiation  into  three  modes.  The  necessary  condition  for  the  oscillator  to  lase  is  that 
the  remaining  terms  exceed  this  loss.  For  the  ultra-high  gain  case  this  requires  ToX  > 
—  ln(-\/ff/3).  This  expression  has  been  confirmed  experimentally13  in  recent  operation 
of  the  NRL  FEL  oscillator.  The  interaction  length,  L,  can  be  varied  by  dumping  the 
beam  at  different  axial  locations  within  the  wiggler.  For  the  following  set  of  experimental 
parameters,  beam  energy  E0  =  500  keV,  beam  current  I  =  100  .4,  wiggler  field  strength 
Bw  =  615  G,  beam  radius  r&  =  0.64  cm  and  wiggier  length  tw  =  4.0  cm,  .he  minimum 
interaction  length  is  determined  to  be  45  cm.  Inserting  this  value  into  Eq.(12)  yields  a 
theoretical  value  of  0.64  for  the  reflection  coefficient.  The  independently  measured  Bragg 
reflection  coefficient  has  the  value  0.65,  which  is  in  excellent  agreement  with  the  theoretical 
value. 
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IV.  Multi-Mode  Simulation 


The  space-time  evolution  of  the  fields  in  the  resonator  is  simulated  by  numerically 
evolving  the  equations  for  the  fields  and  particles.  The  radiation  field  model  for  the  multi- 
mode  simulation  is  given  by, 

4s(z,f)  =  ^an(0  sin  (ft„z)exp(iw»t)el  +  c.c.,  (15) 

n 

<£(M)  =  +k*>)Z  ~  unt\  +  fon{t)  COs[(kn  +  kw)z  -  Unt],  (16) 

n 

where  kn  —  nv/L  =  u>n/c  and  the  sum  is  over  the  discrete  number  of  modes  under  con¬ 
sideration.  This  model  has  the  property  that  the  complex  expansion  coefficients  in  the 
harmonic  analysis,  an(Oi|?in(*)»02>.(f)i  are  only  functions  of  time;  which  results  in  ordi¬ 
nary  differential  equations  for  the  particle  and  field  evolution.  This  model  also  has  the 
attribute  that  the  field  boundary  conditions  at  two  perfectly  reflecting  mirrors  is  automat¬ 
ically  satisfied,  and  we  model  the  resonator  losses  heuristically  by  adding  a  damping  term 
to  the  wave  equation,  [d2  jdz1  -  c~2d2jdt2  -  vc~2d/dt\A[i(z,t)  =  4 ttc~1  J±(z,t),  where 
v  ~  u/Q  and  Q  is  the  quality  factor  of  the  resonator.  The  driving  currents  and  charge  den¬ 
sities  for  the  vector  and  scalar  potentials  are  modeled  with  a  discrete  distribution  function 
as  follows, 


p(z,t)=-e  J  dz0no(t)5{z  -  z(z0,t)),  (17) 

J±{z,  t)  =  ~~  f  dzono(t)5{z  -  z{z0,  t))/~>0,  (18) 

where,  A  =  Ar  +  Aw,  n0(t)  =  n0(cc)(l  -exp(-f/ffl)]  and  n0(oc)  is  the  flattop  density 
of  the  electron  beam  pulse,  Ir  is  the  characteristic  rise  time  for  the  beam  current  or  density, 
and  r(roi  f)  is  the  axial  orbit  of  a  particle  located  at  position  z0  at  t  =  0. 

The  slowly  varying  field  approximation,  dan{t)/dt  <<  yna„(i),  yields  the  following 
set  of  equations  for  the  evolution  of  the  fields  and  particles. 
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A'n  -  3m  f  drQ  sin[ft  (fo,  r)j  P(r0), 

2  -Yo 

(19) 

=  ~~r  +  “p  f  dr0cos[^„(ro,r)]^P(r0), 

2  An  Jt-i  7o 

(20) 

4>m  =  -fon[  dr0cos[ft(r0,r)]P(ro), 

(21) 

Jr-l 

*  r 

fan  = -fan  drosin[0n(ro,r)jP(ro), 

Jt-i 

(22) 

ft  =  +Mi[^mc°s[^m(ro,r)l-^m  sin[^m (r0,  r)J] 

m 

+  P*n  H  {  [(**»  +  k™)L  ~  kmL&:  +  £oz/Mm]Am  sin[0m(ro,  r)j 

-  fozPz&m  cos[ttm(r0,  r)]|.  (23) 

We  have  introduced  the  following  normalized  parameters,  r  =  vQztjL,  An  =  eAn/mc2, 
<p  =  e^/mc2,  (...)'  —  d(...)/dr.  We  have  also  made  the  following  definitions, 

dn  =  — i.4„(r)exp(i0nr),  ft (r0,r)  =  (&n  +  Aw)r(r0,  r)  -  wnXr/v0z  -  0n(r),  = 

uL/voz,  vj  =  w}*.IF(l  -  3l/2)j2klcH  3m  =  F3M*u2p)2c*  3$zkn,  3-m  = 
2Fuj3o/c2(kn  -i-  K)‘13oz,  3m  =  X(*n  +  *w)/T7«^x,  /^n  =  An (kn  +  kw)Lj/23oz7o 

and  P(r0)  =  1  -  exp(-r0/rR),  where  L  is  the  length  of  the  resonator,  F  is  the  filling  factor 
and  wp  is  the  nonrelativist ic  plasma  frequency. 

This  system  of  equations  is  solved  numerically  by  using  a  four-point  Adams-Bashforth 
predictor  corrector  scheme  which  is  initialized  by  using  the  three  point  Runge-Kutta 
method.  The  ensemble  average  over  initial  electrons,  /Tr_1(...)dro,  is  typically  performed 
with  two  thousand  (2000)  particles.  The  results  of  the  simulations  and  the  linear  theory, 
obtained  from  the  linearization  of  Eqs.  (19)  -  (23),  are  shown  in  Figs.  2  and  3. 

V.  Conclusions 

A  comparison  of  the  temporal  growth  rates  obtained  from  the  linear  theory  and  the 
numerical  simulations  is  shown  in  Figs.  2  and  3.  The  growth  rates  for  the  simulations  are 
obtained  numerically  from  the  field  amplitude  data  during  the  initial  field  evolution,  where 
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the  wave  growth  is  linear.  In  Fig.  2  the  data  is  presented  for  a  low  gain  case  with  physical 
parameters  given  by,  7  =  2.0,/  =  5 A,F  =  0.2,  kwrt,  =  0.62831,  ^  =  0.2  a JidLJtw  =  50. 
There  is  an  excellent  agreement  between  theory  and  simulation.  In  Fig.  2  we  also  compare 
the  theoretical  and  numerical  efficiencies  for  a  high  gain  case.  Where  the  efficiency  is 
defined  as  the  stored  electromagnetic  energy  density  normalized  to  the  incident  beam 
energy  density.  The  theoretical  estimates  of  the  efficiency  are  based  on  particle  trapping 
arguments12,  with  the  assumption  that  all  the  energy  lost  by  the  particles  is  converted 
into  electromagnetic  energy.  The  characteristic  change  in  velocity  of  a  particle  is  given  by 
the  difference  in  the  beam  velocity  and  the  phase  velocity  of  the  trapping  potential  This 
phase  velocity  is  approximated  from  the  results  of  the  linear  dispersion  relation.  Again  we 
find  good  qualitative  and  quantitative  agreement  between  the  simulation  and  theory. 
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Comparison  of  the  temporal  growth  rate  from  theory  and  numerical  simulation 
iu  Compton  regime.  The  following  physical  parameters  are  used:  -y  =  2.0, 
"  0.2,  kwrb  =  0.62831,  $w  -  0.2  and  L/tw  =  50. 
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Appendix  A 


In  the  following,  we  shall  consider  the  space-time  evolution  of  the  radiation  fields 
produced  by  the  interaction  of  a  beam  of  relativistic  electrons  with  a  helical  wiggler  field 
contained  within  the  mirrors  of  an  optical  resonator.  The  analysis  is  fully  relativistic 
and  is  conducted  self-consistently  within  the  framework  of  the  Vlasov-Maxwell  system  of 
equations. 

The  wiggler  vector  potential  is  modeled  as  follows, 


B 


Av(z)  =  [exp(i/rw  z)e_  +  exp(-«Aws)e+], 


(.41) 


where  the  wiggler  magnetic  field  strength  is  the  wiggler  period  is  tw  —  jk w  and 
the  basis  vectors  are  e±  =  (e3  ±  ie y)/2.  We  have  assumed  in  this  model  that  the  beam 
radius  is  small  compared  to  the  wiggler  period  (k^ri,  <  l)  hence  the  transverse  gradients 
in  the  wiggler  field  are  neglected.  We  similarly  invoke  the  para-axial  approximation  to  the 
radiation  fields  and  neglect  transverse  coordinate  dependencies  in  the  fields,  to  obtain  the 
following  radiation  field  model  for  the  vector  and  scalar  potentials, 


A{z,t)  =  [d/(z,  t)exp(-ik0z)  +  ab{z,t)exp(ik0z)]exp(iu0t)e-  +  c.c.,  (.42) 

<p{z,  f)  =  o(z,  f)exp[-i(k0  +  kw)z  +  iu0l\  +  c.c.,  (A3) 

where  a/(z, /)  and  a&(z,  t)  denote  the  forward  and  backward  components  of  the  wave  field 
respectively,  and  Uq  =  cko  is  the  frequency.  These  field  coefficients  are  assumed  to  be 
slowly  varying  functions  of  space  and  time  compared  to  the  radiation  wavelength  and 
temporal  period.  The  slow  spatial  dependence  of  the  field  coefficients  is  expressed  by. 
|  Q~ldQ/dz  |<<  ko,  with  Q  =  a/(z, f),  ab(z,  <),o(z,  t)  and  the  slow  temporal  dependence 
of  the  the  field  coefficients  is  expressed  by,  |  Q~ldQ/dt  |<  <  w0- 

The  space-time  evolution  of  the  fields  is  governed  by  Maxwell’s  equations,  which  can 
be  cast  in  the  form,  [d2[dz2  -  c~2d2  jdt2)A[z,t)  =  izc~l  J±{z,  t)  and  d2 /d z2<?(z,  t)  = 
—4 zp(z,t).  The  driving  current  and  charge  densities  are  obtained  from  the  appropriate 
moments  of  the  Vlasov  distribution  function.  The  Vlasov  distribution  function  is  evolved 
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according  to  the  equation,  {djdt  +  [Pz/m^d jdz  -  e[E  +  (px  B)/m7c]  •  d/dp}g(z,p,  t)  = 
0.  By  making  use  of  the  fact  that  the  canonical  transverse  momentum  is  an  invari¬ 
ant  of  the  motion  and  assuming  that  the  beam  is  cold  in  the  transverse  direction  (e.g., 
g{z,Px,Py,pz,t)  =  g(z,pz,  t)8(Px)f>[Py)  )  the  evolution  of  the  reduced  distribution  func¬ 
tion  is  governed  by, 


{ 


Pz  d 
itvjt  ds 


+ 


e2  d 
2m7rc2  dz 


g{z,Pz-.t)  =  o, 


(.44) 


where  mcr^r  =  { m2c *  +c2p2  +  e?(A-  Since  d(Aw  • Aw)/dz  =  0,  the  equilibrium 

distribution  function  satisfies  the  equation,  {d/dt+(pz/m^o)d /dz}g^{z,pz,  t)  =  0,  where 
me2 7 o  =  {m2c*  +  cPp^  -f  e^B^/k^,}1^2.  For  long  electron  beam  pulses  we  shall  consider 
spatially  and  temporally  homogeneous  equilibria  given  by  g^(z,pz,  t)  =  g^{pz).  To  first 
order  in  the  perturbed  fields,  the  evolution  of  the  linearized  distribution  function  is  given 
by, 


dm  e2  d  -  - 
—  e— — I - —  —  (.4  •  .4W) 


dz  m70c2  dz 

The  solution  to  the  linearized  Vlasov  equation  is  formally  given  by, 


dg(0) 

J  dPz 


(.45) 


+  'Xu'J + 


dgw 

dpz 


(.46) 


Linearizing  the  wave  equations  for  .4  and  <t>  one  obtains, 


g{1)(z,P:J) 


m7o 


dp 


<7(0)(P--)  ,  e2  r 


r  "f"  '4a,  (  „ 

m7o  c  mtc 


->/ 


dpz 


9{0){P:) 

m7o 


(.47) 


=  4-n0e  j  dpzg{1)(:,p;,t). 


(.48) 
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By  making  use  of  the  slowly  varying  coefficient  approximations,  the  components  of 
Maxwell’s  equations  can  be  expressed  as: 


(A9) 

(A10) 


{tj  ~  =  *^jjj^«xp(i(Jifz  ~  «o01  J  dP:9{1)[z,pz,  t),  (All) 

where  mc27  =  {mV  +  <rp20  +  e2^2 /** }i/2,  “(n)  =  f  dpz{l/~!o)n9W{Pz), 
w2  =  45rnoe2/rr,  /?„,  =  eBwjm7fc2k-tt]  and  ft  =  k0  +  Also  note  that  in  the  fre¬ 
quency  regime  for  resonant  interaction  of  the  radiation  field  and  the  beam  particles, 
(u'o  »  2^-7? k^c)  the  driving  current  for  the  backward  wave  is  negligible.  We  shall  solve 
this  system  of  equations  in  the  time  asymptotic  limit  for  which  the  forward  wave  oscillates 
at  a  single  complex  frequency,  «/(r,  t)  =s  a/(z)exp(iAwt),  where  Aw  is  to  be  determined 
self-consistently  from  the  boundary  conditions  at  the  mirror  surfaces. 

The  ponderomotive  potential  in  terms  of  the  time  asymptotic  field  coefficients  is  given 
by, 


Aw  ■  A  =  af(z)zxp{-i[Kz  -  (w0  +  Aw)t]}  +  nonresonant  terms.  (A12) 

Retaining  only  the  phase  resonant  terms,  the  formal  solution  to  the  linearized  Vlasov 
equation  can  be  written  as  follows, 


»«(.-, p..,o  =  /*f=a«p{  -UK- 

Jo  Pz  Vz  dpz 


(d/dz’  -  iK)  [^r-aj(z’)  -  0{z’)] . 
-7o 


(-413). 


Inserting  this  result  into  the  Vlasov-Maxwell  system  of  equations  one  obtains 


.  Up  ~  £jo)  (^’O  +  fZ 


2w0<n  &L 


f  dz'(z'  -  ar)exp{  -  i£K(z'  -  r)} 

Jo 


ws.iKn’fwm+Sik 

f*  dz'exp{  -  i&K{z'  -  z)}(d/dz '  -  »A')[(1  +  0%)0v,a/(z')f2  -  fa)],  (.414) 
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where  we  have  defined  ;\K  —  K  —  (u/q  +  &»j)(v:q.  The  previous  set  of  equations  yields  a 
dispersion  relation  which  we  shall  refer  to  as  the  complete  dispersion  relation.  A  simplified 
dispersion  relation  is  obtained  by  noting  that  in  the  momentum  integration  ,  the  results 
are  most  sensitive  to  changes  in  the  exponent,  AK(z'  -  z).  Retaining  only  the  terms  in 
the  integration  by  parts  which  are  proportional  to  d±Kjdpz,  one  obtains  the  following 
simplified  system  of  equations, 
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In  both  the  cases  of  the  complete  and  simplified  set  of  equations,  the  equations  are  of  the 
convolution  type  and  can  be  solved  by  Laplace  transform  methods.  The  text  of  the  paper 
consists  of  a  detailed  analysis  of  the  simplified  set  of  equations  in  the  Compton  regime. 
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